Merging RNAseq with ChIP-seqΒΆ
Dataset Generation for AP-axis analysis¶
What are the genetic and regulatory mechanisms that drive the intense proliferation from the spinal-cord to the forebrain in bilateria?
We know that the silencing of the Hox genes in the forebrain facilitates the growth of the tissue, however, do the Hox genes act in isolation? Or are there other genes that exhibit the same phenotype and thus, may be equally integral to this process?
By removing the silencing mechanism (H3K27me3) we investigate which genes, act in a similar manner by integrating the RNAseq response to the knockout with epigenetic data.
Cells appendix:¶
1) Imports
2) RNAseq merging files
3) Add annotation information
4) Save dataframe combinations to csv files for DE analysis
5) Run R scripts a. DeSEQ2 & b. Normalisation (in R) plot results from normalisation
6) Plot results from DEseq2
7) Make bar chart of all comparisons (i.e. how many DE genes in each comaprison)
8) Select log2(TMM + 1) as normalisation and build merged dataframe (raw + DEseq2 results)
9) Visualise the merged data top most significant genes (Heatmap and PCA)
10) Add histone modification data from Encode (downloaded & processed in APaxis_encodeDataDownload)
11) Filter the merged dataframe (i.e. for significant results)1) Imports¶
In [1]:
"""
--------------------------------------------------------
Import RNAseq data & merge featureCounts files
--------------------------------------------------------
"""
import os, sys
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
import venn
from matplotlib_venn import venn3
from sciviso import *
from sciutil import SciUtil
u = SciUtil()
module_path = os.path.abspath(os.path.join('..'))
sys.path.append(module_path)
"""
--------------------------------------------------------
Global variables
--------------------------------------------------------
"""
date = '20210124'
gene_id = 'entrezgene_id'
gene_name = 'external_gene_name'
sci_colour = ['#483873', '#1BD8A6', '#B117B7', '#AAC7E2', '#FFC107', '#016957', '#9785C0',
'#D09139', '#338A03', '#FF69A1', '#5930B1', '#FFE884', '#35B567', '#1E88E5',
'#ACAD60', '#A2FFB4', '#B618F5', '#854A9C']
hist_colour = '#483873'
e11_colour = 'white'
e13_colour = 'lightgrey'
e15_colour = '#A7A9AC'
e18_colour = '#58595B'
fb_colour = '#ffbf80'
mb_colour = '#ff8c1a'
hb_colour = '#b35900'
sc_colour = '#663300'
grey = '#bdbdbd'
h3k36me3_colour = '#CEF471'
h3k27me3_colour = '#9F71F4'
h3k4me3_colour = '#9F00FA'
h3k4me2_colour = '#5930B1'
h3k4me1_colour = '#FFE884'
h3k27ac_colour = '#35B567'
h3k9me3_colour = '#1E88E5'
h3k9ac_colour = '#A2FFB4'
wt_colour = '#AADFF1'
ko_colour = '#A53736'
sns.palplot(sci_colour)
sns.color_palette(sci_colour)
project_name = 'prelim'
data_dir = '../data/'
r_dir = f'{data_dir}results/deseq2/'
fig_dir = f'../figures/{project_name}/'
output_dir = f'{data_dir}results/{project_name}/'
input_dir = f'{data_dir}input/'
rna_dir = f'{input_dir}feature_counts/'
cmap = 'RdBu_r'
u_id = 'u_id'
sns.palplot(sci_colour)
sns.color_palette(sci_colour)
def get_time_colour(c):
if '11' in c or '10' in c:
return e11_colour
elif '13' in c or '12' in c:
return e13_colour
elif '15' in c or '14' in c:
return e15_colour
elif '18' in c or '16' in c:
return e18_colour
return '#FFFFFF'
def get_tissue_colour(c):
if 'sc' in c or 'spinal' in c:
return sc_colour
elif 'hb' in c or 'hindbrain' in c:
return hb_colour
elif 'mb' in c or 'midbrain' in c:
return mb_colour
elif 'fb' in c or 'forebrain' in c:
return fb_colour
return '#FFFFFF'
def get_cond_colour(c):
if 'ko' in c:
return ko_colour
elif 'wt' in c:
return wt_colour
return '#FFFFFF'
def get_mark_colour(c):
if '36me3' in c:
return h3k36me3_colour
elif '27me3' in c:
return h3k27me3_colour
elif 'K4me3' in c:
return h3k4me3_colour
elif 'K4me2' in c:
return h3k4me2_colour
elif 'K4me1' in c:
return h3k4me1_colour
elif '27ac' in c:
return h3k27ac_colour
elif 'K9me3' in c:
return h3k9me3_colour
elif 'K9ac' in c:
return h3k9ac_colour
return '#FFFFFF'
def pplot():
plt.rcParams['figure.figsize'] = [5, 4]
plt.rcParams['image.cmap'] = cmap
plt.rcParams['svg.fonttype'] = 'none'
def cplot():
plt.rcParams['figure.figsize'] = [6, 4]
sns.color_palette(sci_colour)
def save_fig(title, ending='.svg'):
plt.savefig(f'{fig_dir}{title.replace(" ", "-")}{ending}')
2) RNAseq merging files¶
Here we import the results from our FeatureCounts data. We merge the files together on their entrez ID.
In [2]:
files = os.listdir(rna_dir)
rna_files = []
for f in files:
if 'summary' not in f:
rna_files.append(f)
# Make a df out of all expression data
rna_files.sort()
df = pd.DataFrame()
# Basically want to create a dummy df since we don't want the columns except the expression
# and the entrez gene id
tmp_df = pd.read_csv(f'{rna_dir}{rna_files[0]}', sep='\t', header=1)
df['Geneid'] = tmp_df['Geneid'].values
df[rna_files[0][:-4]] = tmp_df[tmp_df.columns[-1]].values
# Now we'll add all the other RNA files to this dataframe by joining on ID
for filename in rna_files[1:]:
u.dp(["Adding", filename])
tmp_df = pd.DataFrame()
file_df = pd.read_csv(f'{rna_dir}{filename}', sep='\t', header=1)
print(file_df.head())
u.dp([f'Length of {filename}: ', len(file_df)])
tmp_df['Geneid'] = file_df['Geneid'].values
# Note we take the log transform here
tmp_df[filename[:-4]] = file_df[file_df.columns[-1]].values
df = df.merge(tmp_df, on='Geneid', how='outer')
# Save DF to csv for normalisation
u.dp(["Length of merged RNAseq dataframe: ", len(df)])
3) Add annotation information¶
We add annotation information to the genes: name, etc. This will also allow us to merge with our ChIP data which have been annotated with ensembl IDs.
We also want to save the merged dataframe to CSV files (grouped by the tissue) as then we will run DeSEQ2 on these count files.
In [3]:
"""
--------------------------------------------------------
Add annotation to RNAseq dataframe & save Tissue counts as files.
--------------------------------------------------------
"""
# Lastly let's add some gene annotation information
"""
Annotations were generated using scibiomart v1.0.0 CLI
scibiomart --m ENSEMBL_MART_ENSEMBL --d mmusculus_gene_ensembl --a "ensembl_gene_id,entrezgene_id,external_gene_name,chromosome_name,start_position,end_position,strand" --o mm10Sorted_ --s t
"""
# Remove data that were not able to be annotated to gene names (i.e. these are usually the genes which are "unknown")
annot_df = pd.read_csv(os.path.join(input_dir, 'supps', 'mm10Sorted_mmusculus_gene_ensembl-GRCm38.p6.csv'))
gene_names = annot_df[gene_name].values
gene_id_to_name = {}
for i, g in enumerate(annot_df[gene_id].values):
if not gene_id_to_name.get(g):
gene_id_to_name[g] = gene_names[i]
# add these to the df column
gene_names = []
pseudo_id = []
df = df.rename(columns={'Geneid': 'entrezgene_id'})
for g in df[gene_id].values:
gene_n = gene_id_to_name.get(g)
gene_names.append(gene_n)
pseudo_id.append(f'{g}-{gene_n}')
df[gene_name] = gene_names
df['u_id'] = pseudo_id
# Organise columns a bit nicer
cols = [gene_id, gene_name, 'u_id',
'ko11fb1', 'ko11fb2', 'ko13fb1', 'ko13fb2', 'ko15fb1', 'ko15fb2', 'ko18fb1', 'ko18fb2',
'ko11mb1', 'ko11mb2', 'ko13mb1', 'ko13mb2', 'ko15mb1', 'ko15mb2', 'ko18mb1', 'ko18mb2',
'ko11hb1', 'ko11hb2', 'ko13hb1', 'ko13hb2', 'ko15hb1', 'ko15hb2', 'ko18hb1', 'ko18hb2',
'ko11sc1', 'ko11sc2', 'ko13sc1', 'ko13sc2', 'ko15sc1', 'ko15sc2', 'ko18sc1', 'ko18sc2',
'wt11fb1', 'wt11fb2', 'wt13fb1', 'wt13fb2', 'wt15fb1', 'wt15fb2', 'wt18fb1', 'wt18fb2',
'wt11mb1', 'wt11mb2', 'wt13mb1', 'wt13mb2', 'wt15mb1', 'wt15mb2', 'wt18mb1', 'wt18mb2',
'wt11hb1', 'wt11hb2', 'wt13hb1', 'wt13hb2', 'wt15hb1', 'wt15hb2', 'wt18hb1', 'wt18hb2',
'wt11sc1', 'wt11sc2', 'wt13sc1', 'wt13sc2', 'wt15sc1', 'wt15sc2', 'wt18sc1', 'wt18sc2'
]
df = df[cols]
u.dp(["Number of genes in total:", len(df)])
# Remove the values that don't map to a gene name. On investigation these are
# most commonly predicted genes, or non-coding RNA (which we don't expect to have messenger RNA for)
# for example: 2610203C22Rik RIKEN cDNA 2610203C22 gene [ Mus musculus (house mouse) ]
# There were about 6000 genes. None of which had a logFC > 2 with the majority having 0 LogFC
df = df[df[gene_name].values != None]
u.dp(["Number of genes with gene names:", len(df)])
df.to_csv(f'{output_dir}merged_df_FEATURE_COUNTS_annot_{date}.csv', index=False)
4) Save dataframe combinations to csv files¶
Here we just save each combination that we want to run DE on to a CSV file.
In [5]:
"""
Save dataframe combinations to CSV files to do DEseq2 analysis
"""
# Create DF for each tissue
tissues = ['fb', 'mb', 'hb', 'sc']
for t in tissues:
ts = [u_id]
for c in df.columns:
if t in c and '11' not in c and 'id' not in c:
ts.append(c)
t_df = df[ts]
t_df.to_csv(f'{r_dir}merged_df_{t}_FEATURE_COUNTS_{date}.csv', index=False)
# Create DF for each timepoint
times = ['11', '13', '15', '18']
for t in times:
ts = [u_id]
for c in df.columns:
if t in c and 'id' not in c and ('fb' in c or 'mb' in c):
ts.append(c)
t_df = df[ts]
t_df.to_csv(f'{r_dir}merged_df_anterior_{t}_FEATURE_COUNTS_{date}.csv', index=False)
for t in times:
ts = [u_id]
for c in df.columns:
if t in c and 'id' not in c and ('sc' in c or 'hb' in c):
ts.append(c)
t_df = df[ts]
t_df.to_csv(f'{r_dir}merged_df_posterior_{t}_FEATURE_COUNTS_{date}.csv', index=False)
# Create df for all the data (including e11)
columns = [gene_id]
for c in df.columns:
if 'wt' in c or 'ko' in c:
columns.append(c)
df[columns].to_csv(f'{r_dir}merged_df_FEATURE_COUNTS_{date}.csv', index=False)
# Create combo between each tissue (i.e sc vs mb)
def gen_csv_combo(cond1, cond2):
condition = ['wt', 'ko']
for t in condition:
ts = [u_id]
for c in df.columns:
if t in c and (cond1 in c or cond2 in c) and gene_id != c and gene_name != c and '11' not in c:
ts.append(c)
t_df = df[ts]
t_df.to_csv(f'{r_dir}merged_df_{t}_{cond1}-{cond2}_FEATURE_COUNTS_{date}.csv', index=False)
u.dp(["Cond done:", f'{t}_{cond1}-{cond2}'])
def gen_csv_combo_time(cond1, cond2):
condition = ['wt', 'ko']
for t in condition:
ts = [u_id]
for c in df.columns:
if t in c and (cond1 in c or cond2 in c) and gene_id != c and gene_name != c and ('sc' in c or 'hb' in c):
ts.append(c)
t_df = df[ts]
t_df.to_csv(f'{r_dir}merged_df_posterior_{t}_{cond1}-{cond2}_FEATURE_COUNTS_{date}.csv', index=False)
u.dp(["Cond done:", f'posterior_{t}_{cond1}-{cond2}'])
for t in condition:
ts = [u_id]
for c in df.columns:
if t in c and (cond1 in c or cond2 in c) and gene_id != c and gene_name != c and ('fb' in c or 'mb' in c):
ts.append(c)
t_df = df[ts]
t_df.to_csv(f'{r_dir}merged_df_anterior_{t}_{cond1}-{cond2}_FEATURE_COUNTS_{date}.csv', index=False)
u.dp(["Cond done:", f'anterior_{t}_{cond1}-{cond2}'])
gen_csv_combo('fb', 'mb')
gen_csv_combo('fb', 'hb')
gen_csv_combo('fb', 'sc')
gen_csv_combo('mb', 'sc')
gen_csv_combo('mb', 'hb')
gen_csv_combo('hb', 'sc')
gen_csv_combo_time('11', '18')
gen_csv_combo_time('11', '13')
gen_csv_combo_time('13', '18')
gen_csv_combo_time('15', '18')
gen_csv_combo_time('11', '15')
gen_csv_combo_time('13', '15')
df.to_csv(f'{r_dir}merged_df_FEATURE_COUNTS_annot_{date}.csv', index=False)
u.dp(["Combinations generated!"])
5) Run R scripts a. DeSEQ2 & b. Normalisation¶
Here we need to step out of this notebook and run the following R markdowns:
2) Run Deseq2 (apAxis_between-cond-tissue.Rmd, apAxis_between-time.Rmd, apAxis_between-tissue.Rmd)
3) Normalise the data (apAxis_between-cond-tissue.Rmd)
Prior to running our analysis, we need to normalise the RNAseq data. Given we're interested in integrating the changes across the genes, we use EdgeRs TMM method.
Script below for normalisation¶
Run script in file: apAxis_between-cond-tissue.RmdIn [6]:
"""
--------------------------------------------------------
Display the normalisations
--------------------------------------------------------
"""
%matplotlib inline
from sciviso import Histogram, Scatterplot
plot_on = True
plt_on = True
save_on = True
# Show an example of each method's normalisation
def show_scatter(df, cond_x, cond_y, x_label='', y_label='', title='', log2=False, offset=1, annotations=None):
a = 0.1
offset = 1
plt_df = pd.DataFrame()
x = cond_x
y = cond_y
if log2:
x = f'Log2({cond_x} + {offset})'
x = f'Log2({cond_y} + {offset})'
plt_df[x] = np.log2(df[cond_x].values + offset)
plt_df[y] = np.log2(df[cond_y].values + offset)
else:
plt_df[x] = np.log2(df[cond_x].values + offset)
plt_df[y] = np.log2(df[cond_y].values + offset)
# plot histogram
h = Histogram(plt_df, x, title=title, fit_norm=False, plot_rug=True, xlabel=x)
h.plot()
save_fig(f'Hist_{x}_{title}', ending='.pdf')
plt.show()
# now run the scatter
s = Scatterplot(plt_df, x, y, title=title, xlabel=x, ylabel=y, colour="grey", add_legend=False)
s.opacity=0.5
s.plot()
save_fig(f'Scatter_{x}-{y}_{title}', ending='.pdf')
plt.show()
deseq2 = pd.read_csv(f'{r_dir}merged_df_FEATURE_COUNTS_DEseq2Norm_{date}.csv')
rlog = pd.read_csv(f'{r_dir}merged_df_FEATURE_COUNTS_rlog_{date}.csv')
tmm = pd.read_csv(f'{r_dir}merged_df_FEATURE_COUNTS_tmm_{date}.csv')
vst = pd.read_csv(f'{r_dir}merged_df_FEATURE_COUNTS_vst_{date}.csv')
sc_genes = None
show_scatter(deseq2, 'wt11fb1', 'ko11fb1', 'wt11fb1', 'ko11fb1', 'Forebrain 11 WT vs KO DEseq2', True, 1.0, sc_genes)
show_scatter(deseq2, 'wt11fb1', 'wt11fb2', 'wt11fb1', 'wt11fb2', 'Forebrain WT vs WT DEseq2', True, 1.0, sc_genes)
show_scatter(rlog, 'wt11fb1', 'ko11fb1', 'wt11fb1', 'ko11fb1', 'Forebrain 11 WT vs KO Rlog', False, 1.0, sc_genes)
show_scatter(rlog, 'wt11fb1', 'wt11fb2', 'wt11fb1', 'wt11fb2', 'Forebrain WT vs WT Rlog', False, 1.0, sc_genes)
show_scatter(tmm, 'wt11fb1', 'ko11fb1', 'wt11fb1', 'ko11fb1', 'Forebrain 11 WT vs KO TMM', True, 1.0, sc_genes)
show_scatter(tmm, 'wt11fb1', 'wt11fb2', 'wt11fb1', 'wt11fb2', 'Forebrain WT vs WT TMM', True, 1.0, sc_genes)
show_scatter(vst, 'wt11fb1', 'ko11fb1', 'wt11fb1', 'ko11fb1', 'Forebrain 11 WT vs KO VST', False, 1.0, sc_genes)
show_scatter(vst, 'wt11fb1', 'wt11fb2', 'wt11fb1', 'wt11fb2', 'Forebrain WT vs WT VST', False, 1.0, sc_genes)
6) Plot results from DEseq2¶
Here we want to do the normal figures for a DE analysis i.e. volcano plots, heatmaps, & Venn diagrams etc.
In [7]:
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from adjustText import adjust_text
from sciviso import Vis
class Volcanoplot(Vis):
def __init__(self, df: pd.DataFrame, log_fc: str, p_val: str, label_column: str, title='',
xlabel='', ylabel='', invert=False, p_val_cutoff=0.05,
log_fc_cuttoff=2, label_big_sig=False, colours=None, offset=None,
text_colours=None, values_to_label=None, max_labels=20, values_colours=None,
figsize=(2, 2), title_font_size=8, label_font_size=6, title_font_weight=700):
super().__init__(df, figsize=figsize, title_font_size=title_font_size, label_font_size=label_font_size,
title_font_weight=title_font_weight)
super().__init__(df)
self.log_fc = log_fc
self.p_val = p_val
self.p_val_cutoff = p_val_cutoff
self.log_fc_cuttoff = log_fc_cuttoff
self.values_to_label = values_to_label
self.label_big_sig = label_big_sig
self.invert = invert
self.label_column = label_column
self.offset = offset
self.label = 'volcanoplot'
self.colours = {'ns_small-neg-logFC': 'lightgrey',
'ns_small-pos-logFC': 'lightgrey',
'ns_big-neg-logFC': 'grey',
'ns_big-pos-logFC': 'grey',
'sig_small-neg-logFC': '#2970b1',
'sig_small-pos-logFC': '#d6604c',
'sig_big-neg-logFC': '#0a3568',
'sig_big-pos-logFC': '#6f0220'} if colours is None else colours
self.xlabel = xlabel
self.ylabel = ylabel
self.title = title
self.max_labels = max_labels
self.values_colours = values_colours or {}
self.text_colours = text_colours or {}
def add_scatter_and_annotate(self, fig: plt, x_all: np.array, y_all: np.array,
colour: str, idxs: np.array, annotate=False):
x = x_all[idxs]
y = y_all[idxs]
ax = fig.scatter(x, y, c=colour, alpha=self.opacity, s=10, vmin=-10.0, vmax=10.0)
# Check if we want to annotate any of these with their gene IDs
if self.values_to_label is not None:
texts = []
labels = self.df[self.label_column].values[idxs]
for i, name in enumerate(labels):
if name in self.values_to_label:
lbl_bg = self.values_colours.get(name)
color = self.text_colours.get(name)
texts.append(fig.text(x[i], y[i], name, color=color, fontsize=6,
bbox=dict(fc=lbl_bg, alpha=1.0, boxstyle='round,pad=0.1', lw=0)))
adjust_text(texts, force_text=2.0)
# Check if the user wants these labeled
if self.label_big_sig and annotate:
# If they do have a limit on the number of ones we show (i.e. we don't want 10000 gene names...)
max_values = -1 * self.max_labels
if len(y) < self.max_labels:
max_values = -1 * (len(y) - 1)
most_sig_idxs = np.argpartition(y, max_values)[max_values:]
labels = self.df[self.label_column].values[idxs][most_sig_idxs]
x = x[most_sig_idxs]
y = y[most_sig_idxs]
# We only label the ones with the max log fc
for i, name in enumerate(labels):
fig.annotate(name, (x[i], y[i]),
xytext=(0, 10),
textcoords='offset points', ha='center', va='bottom',
bbox=dict(boxstyle='round,pad=0',
fc=None, alpha=0.2)
)
return ax
def plot(self):
"""
For annotation styling see: https://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.annotate
Returns
-------
"""
# if offset is not given, make the offset the smallest value in the dataset
if not self.offset:
vals = self.df[self.p_val].values
self.offset = np.min(vals[np.nonzero(vals)])
self.u.warn_p(['No offset was provided, setting offset to be smallest value recorded in dataset: ',
self.offset])
# x axis has log_fc, first only plot the values < cutoff
x = self.df[self.log_fc].values
y = -1 * np.log10(self.df[self.p_val].values + self.offset)
log_fc_np = self.df[self.log_fc].values
p_val_np = self.df[self.p_val].values
if self.invert:
x = -1 * np.log10(self.df[self.p_val].values + self.offset)
y = self.df[self.log_fc].values
sig_small_pos_logfc = np.where((p_val_np <= self.p_val_cutoff) & (np.abs(log_fc_np) < self.log_fc_cuttoff)
& (log_fc_np > 0))
sig_big_pos_logfc = np.where((p_val_np <= self.p_val_cutoff) & (np.abs(log_fc_np) >= self.log_fc_cuttoff)
& (log_fc_np > 0))
sig_small_neg_logfc = np.where((p_val_np <= self.p_val_cutoff) & (np.abs(log_fc_np) < self.log_fc_cuttoff)
& (log_fc_np <= 0))
sig_big_neg_logfc = np.where((p_val_np <= self.p_val_cutoff) & (np.abs(log_fc_np) >= self.log_fc_cuttoff)
& (log_fc_np <= 0))
# Plot the points
fig, ax = plt.subplots(figsize=(2.5, 2.5))
self.add_scatter_and_annotate(ax, x, y, self.colours['sig_small-pos-logFC'], sig_small_pos_logfc)
self.add_scatter_and_annotate(ax, x, y, self.colours['sig_big-pos-logFC'], sig_big_pos_logfc, annotate=True)
# Negative
self.add_scatter_and_annotate(ax, x, y, self.colours['sig_small-neg-logFC'], sig_small_neg_logfc)
self.add_scatter_and_annotate(ax, x, y, self.colours['sig_big-neg-logFC'], sig_big_neg_logfc, annotate=True)
self.add_labels()
ax.tick_params(labelsize=6)
ax.tick_params(direction='out', length=2, width=0.5)
ax.spines['bottom'].set_linewidth(0.5)
ax.spines['top'].set_linewidth(0)
ax.spines['left'].set_linewidth(0.5)
ax.spines['right'].set_linewidth(0)
ax.tick_params(labelsize=6)
ax.tick_params(axis='x', which='major', pad=0)
ax.tick_params(axis='y', which='major', pad=0)
return ax
In [11]:
"""
---------------------------------------------------------------
Read in results from DEseq2 and format DFs nicer
---------------------------------------------------------------
"""
df_dict = {}
deseq2_files = os.listdir(r_dir)
for f in deseq2_files:
if 'DEseq2' in f:
try:
de_df = pd.read_csv(os.path.join(r_dir, f))
de_df = de_df.rename(columns={de_df.columns[0]: 'u_id'})
gene_names = [s.split('-')[1] for s in list(de_df['u_id'].values)]
gene_ids = [s.split('-')[0] for s in list(de_df['u_id'].values)]
de_df['padj'] = de_df['padj'].fillna(1) # Replace Nan p values with 1.0s
# Replace Nans with 0's for other values
de_df = de_df.replace(np.nan, 0)
de_df[gene_id] = gene_ids
de_df[gene_name] = gene_names
df_dict[f] = de_df
except:
print(f)
"""
---------------------------------------------------------------
Set up markers to display on Volcano plot
---------------------------------------------------------------
"""
fb_genes = ['Emx1', 'Eomes', 'Tbr1', 'Foxg1', 'Lhx6']
mb_genes = ['En1', 'En2', 'Lmx1a', 'Bhlhe23', 'Sall4']
hb_genes = ['Phox2b', 'Krox20', 'Fev', 'Hoxb1', 'Hoxd3']
sc_genes = ['Hoxd8', 'Hoxd9', 'Hoxd10', 'Hoxd11','Hoxa7', 'Hoxa9', 'Hoxa10',
'Hoxb9', 'Hoxb13', 'Hoxc8', 'Hoxc9', 'Hoxc10', 'Hoxc11', 'Hoxc12', 'Hoxc13']
progenitors = ['Sox2', 'Sox1', 'Sox3', 'Hes1', 'Hes5']
neurons = ['Tubb3', 'Snap25', 'Syt1', 'Slc32a1','Slc17a6']
glia = ['Pdgfra', 'Cspg4', 'Aqp4', 'Egfr', 'Slc6a11']
vals_to_label = fb_genes + mb_genes + hb_genes + sc_genes + progenitors + neurons + glia
lbl_colours = {}
text_colours = {}
for g in fb_genes:
lbl_colours[g] = fb_colour
text_colours[g] = 'black'
for g in mb_genes:
lbl_colours[g] = mb_colour
text_colours[g] = 'black'
for g in hb_genes:
lbl_colours[g] = hb_colour
text_colours[g] = 'white'
for g in sc_genes:
lbl_colours[g] = sc_colour
text_colours[g] = 'white'
for g in progenitors:
lbl_colours[g] = e11_colour
text_colours[g] = 'black'
for g in neurons:
lbl_colours[g] = e13_colour
text_colours[g] = 'black'
for g in glia:
lbl_colours[g] = e15_colour
text_colours[g] = 'black'
def plot_volcano(df, title, save=True, show=True):
volcanoplot = Volcanoplot(df, 'log2FoldChange', 'padj', gene_name,
title, 'Log 2 Fold change', '-log10(p adj)',
p_val_cutoff=0.05,
label_big_sig=False, log_fc_cuttoff=1.5,
values_to_label=vals_to_label, figsize=(2,2),
values_colours=lbl_colours, text_colours=text_colours)
sns.set_style("ticks")
volcanoplot.plot()
if save:
save_fig(f'Volcano{title.replace(" ", "")}', ending='.pdf')
if show:
plt.show()
"""
---------------------------------------------------------------
Plot volcanos
---------------------------------------------------------------
"""
plot_volcano(df_dict[f'DEseq2_CNS_fb_{date}.csv'], "FB EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_mb_{date}.csv'], "MB EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_hb_{date}.csv'], "HB EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_sc_{date}.csv'], "SC EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_anterior_11_{date}.csv'], "Anterior e11.5 EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_anterior_13_{date}.csv'], "Anterior e13.5 EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_anterior_15_{date}.csv'], "Anterior e15.5 EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_anterior_18_{date}.csv'], "Anterior e18.5 EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_posterior_11_{date}.csv'], "Posterior e11.5 EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_posterior_13_{date}.csv'], "Posterior e13.5 EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_posterior_15_{date}.csv'], "Posterior e15.5 EedcKO vs WT")
plot_volcano(df_dict[f'DEseq2_CNS_posterior_18_{date}.csv'], "Posterior e18.5 EedcKO vs WT")
In [12]:
"""
---------------------------------------------------------------
Plot venn diagrams
---------------------------------------------------------------
"""
import venn
from matplotlib.colors import ListedColormap
def plot_venn(gene_sets, labels, title, save=True, show=True, colours=None):
output = f'venn_{title.replace(" ", "-")}'
dset = {}
for i, l in enumerate(labels):
dset[l] = gene_sets[i]
if len(gene_sets) > 3:
cmap = ListedColormap(sns.color_palette(colours))
venn.venn(dset, cmap=cmap, fontsize=6, figsize=(2, 2))
else:
venn3(gene_sets, set_labels=labels)
if save:
save_fig(output)
if show:
plt.show()
fb_df = df_dict[f'DEseq2_CNS_fb_{date}.csv'][df_dict[f'DEseq2_CNS_fb_{date}.csv']['log2FoldChange'] > 0.5]
mb_df = df_dict[f'DEseq2_CNS_mb_{date}.csv'][df_dict[f'DEseq2_CNS_mb_{date}.csv']['log2FoldChange'] > 0.5]
hb_df = df_dict[f'DEseq2_CNS_hb_{date}.csv'][df_dict[f'DEseq2_CNS_hb_{date}.csv']['log2FoldChange'] > 0.5]
sc_df = df_dict[f'DEseq2_CNS_sc_{date}.csv'][df_dict[f'DEseq2_CNS_sc_{date}.csv']['log2FoldChange'] > 0.5]
fb_sig_genes = fb_df[fb_df['padj'] < 0.05][gene_name].values
mb_sig_genes = mb_df[mb_df['padj'] < 0.05][gene_name].values
hb_sig_genes = hb_df[hb_df['padj'] < 0.05][gene_name].values
sc_sig_genes = sc_df[sc_df['padj'] < 0.05][gene_name].values
plot_venn([set(fb_sig_genes),
set(mb_sig_genes), set(hb_sig_genes), set(sc_sig_genes)],
['Forebrain WT vs KO', 'Midbrain WT vs KO', 'Hindbrain WT vs KO', 'Spinalcord WT vs KO'],
'tissue_WTvsKO', colours=[fb_colour, mb_colour, hb_colour, sc_colour])
plot_venn((set(fb_sig_genes), set(mb_sig_genes), set(sc_sig_genes)),
('Forebrain', 'Midbrain', 'Spinal cord'), 'FbMbSc')
"""
---------------------------------------------------------------
Plot Anterior Venns
---------------------------------------------------------------
"""
e11 = df_dict[f'DEseq2_CNS_anterior_11_{date}.csv'][df_dict[f'DEseq2_CNS_anterior_11_{date}.csv']['log2FoldChange'] > 0.5]
e13 = df_dict[f'DEseq2_CNS_anterior_13_{date}.csv'][df_dict[f'DEseq2_CNS_anterior_13_{date}.csv']['log2FoldChange'] > 0.5]
e15 = df_dict[f'DEseq2_CNS_anterior_15_{date}.csv'][df_dict[f'DEseq2_CNS_anterior_15_{date}.csv']['log2FoldChange'] > 0.5]
e18 = df_dict[f'DEseq2_CNS_anterior_18_{date}.csv'][df_dict[f'DEseq2_CNS_anterior_18_{date}.csv']['log2FoldChange'] > 0.5]
e11_sig_genes = e11[e11['padj'] < 0.05][gene_name].values
e13_sig_genes = e13[e13['padj'] < 0.05][gene_name].values
e15_sig_genes = e15[e15['padj'] < 0.05][gene_name].values
e18_sig_genes = e18[e18['padj'] < 0.05][gene_name].values
plot_venn([set(e11_sig_genes),
set(e13_sig_genes), set(e15_sig_genes), set(e18_sig_genes)],
['e11 WT vs KO', 'e13 WT vs KO', 'e15 WT vs KO', 'e18 WT vs KO'],
'posterior_time_WTvsKO', colours=[e11_colour, e13_colour, e15_colour, e18_colour])
"""
---------------------------------------------------------------
Plot posterior venn diagrams
---------------------------------------------------------------
"""
e11 = df_dict[f'DEseq2_CNS_posterior_11_{date}.csv'][df_dict[f'DEseq2_CNS_posterior_11_{date}.csv']['log2FoldChange'] > 0.5]
e13 = df_dict[f'DEseq2_CNS_posterior_13_{date}.csv'][df_dict[f'DEseq2_CNS_posterior_13_{date}.csv']['log2FoldChange'] > 0.5]
e15 = df_dict[f'DEseq2_CNS_posterior_15_{date}.csv'][df_dict[f'DEseq2_CNS_posterior_15_{date}.csv']['log2FoldChange'] > 0.5]
e18 = df_dict[f'DEseq2_CNS_posterior_18_{date}.csv'][df_dict[f'DEseq2_CNS_posterior_18_{date}.csv']['log2FoldChange'] > 0.5]
e11_sig_genes = e11[e11['padj'] < 0.05][gene_name].values
e13_sig_genes = e13[e13['padj'] < 0.05][gene_name].values
e15_sig_genes = e15[e15['padj'] < 0.05][gene_name].values
e18_sig_genes = e18[e18['padj'] < 0.05][gene_name].values
plot_venn([set(e11_sig_genes),
set(e13_sig_genes), set(e15_sig_genes), set(e18_sig_genes)],
['e11 WT vs KO', 'e13 WT vs KO', 'e15 WT vs KO', 'e18 WT vs KO'],
'posterior_time_WTvsKO', colours=[e11_colour, e13_colour, e15_colour, e18_colour])
"""
---------------------------------------------------------------
Print out the number of significant genes
---------------------------------------------------------------
"""
fb_sig_genes = fb_df[fb_df['padj'] < 0.05][gene_name].values
mb_sig_genes = mb_df[mb_df['padj'] < 0.05][gene_name].values
hb_sig_genes = hb_df[hb_df['padj'] < 0.05][gene_name].values
sc_sig_genes = sc_df[sc_df['padj'] < 0.05][gene_name].values
u.dp(["P.adj < 0.05 for SC: KO vs WT", len(sc_sig_genes)])
u.dp(["P.adj < 0.05 for HB: KO vs WT", len(hb_sig_genes)])
u.dp(["P.adj < 0.05 for MB: KO vs WT", len(mb_sig_genes)])
u.dp(["P.adj < 0.05 for FB: KO vs WT", len(fb_sig_genes)])
7) Make bar chart of all comparisons¶
Have a look at the comaprisons and just plot the number of DEGs for a summary figure.
In [13]:
"""
---------------------------------------------------------------
Bar charts for tissue specific DE
---------------------------------------------------------------
"""
dfs = [df_dict[f'DEseq2_CNS_wt_fb-mb_{date}.csv'], df_dict[f'DEseq2_CNS_ko_fb-mb_{date}.csv'],
df_dict[f'DEseq2_CNS_wt_mb-hb_{date}.csv'], df_dict[f'DEseq2_CNS_ko_mb-hb_{date}.csv'],
df_dict[f'DEseq2_CNS_wt_hb-sc_{date}.csv'], df_dict[f'DEseq2_CNS_ko_hb-sc_{date}.csv'],
df_dict[f'DEseq2_CNS_wt_mb-sc_{date}.csv'], df_dict[f'DEseq2_CNS_ko_mb-sc_{date}.csv'],
df_dict[f'DEseq2_CNS_wt_fb-hb_{date}.csv'], df_dict[f'DEseq2_CNS_ko_fb-hb_{date}.csv'],
df_dict[f'DEseq2_CNS_wt_fb-sc_{date}.csv'], df_dict[f'DEseq2_CNS_ko_fb-sc_{date}.csv'],
df_dict[f'DEseq2_CNS_fb_{date}.csv'],
df_dict[f'DEseq2_CNS_mb_{date}.csv'],
df_dict[f'DEseq2_CNS_hb_{date}.csv'],
df_dict[f'DEseq2_CNS_sc_{date}.csv']
]
labels = ['Forebrain vs Midbrain WT',
'Forebrain vs Midbrain KO',
'Midbrain vs Hindbrain WT',
'Midbrain vs Hindbrain KO',
'Hindbrain vs Spinalcord WT',
'Hindbrain vs Spinalcord KO',
'Midbrain vs Spinalcord WT',
'Midbrain vs Spinalcord KO',
'Forebrain vs Hindbrain WT',
'Forebrain vs Hindbrain KO',
'Forebrain vs Spinal cord WT',
'Forebrain vs Spinalcord KO',
'WT vs KO (Forebrain)',
'WT vs KO (Midbrain)',
'WT vs KO (Hindbrain)',
'WT vs KO (Spinal cord)',
]
i = 0
values = []
logfc_cutoff = 0.5
for d in dfs:
sigD = d[d['padj'] < 0.05]
sigD = sigD[abs(sigD['log2FoldChange']) > logfc_cutoff]
num_sig = len(sigD)
values.append(num_sig)
print(f'{num_sig}, {labels[i]}')
i += 1
# Plot a simple bar chart.
barchart = Barchart(df, labels, values, order=labels, figsize=(1,1))
barchart.palette=sns.color_palette([wt_colour, ko_colour, wt_colour, ko_colour, wt_colour, ko_colour,
wt_colour, ko_colour, wt_colour, ko_colour, wt_colour, ko_colour, fb_colour, mb_colour, hb_colour, sc_colour])
barchart.figsize=(2,1)
barchart.plot()
save_fig(f'Barchart_numSigs')
u.dp(["Barchart used p.adj cutoff:", 0.05, "\nLogFC cutoff:", logfc_cutoff])
In [14]:
"""
---------------------------------------------------------------
Bar charts for time specific DE (anterior)
---------------------------------------------------------------
"""
dfs = [df_dict[f'DEseq2_CNS_anterior_wt_11-13_{date}.csv'], df_dict[f'DEseq2_CNS_anterior_ko_11-13_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_wt_13-15_{date}.csv'], df_dict[f'DEseq2_CNS_anterior_ko_13-15_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_wt_15-18_{date}.csv'], df_dict[f'DEseq2_CNS_anterior_ko_15-18_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_wt_11-15_{date}.csv'], df_dict[f'DEseq2_CNS_anterior_ko_11-15_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_wt_13-18_{date}.csv'], df_dict[f'DEseq2_CNS_anterior_ko_13-18_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_wt_11-18_{date}.csv'], df_dict[f'DEseq2_CNS_anterior_ko_11-18_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_11_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_13_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_15_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_18_{date}.csv']
]
labels = ['e11 vs e13 (WT)', 'e11 vs e13 (KO)',
'e13 vs e15 (WT)', 'e13 vs e15 (KO)',
'e15 vs e18 (WT)', 'e15 vs e18 (KO)',
'e11 vs e15 (WT)', 'e11 vs e15 (KO)',
'e13 vs e18 (WT)', 'e13 vs e18 (KO)',
'e11 vs e18 (WT)', 'e11 vs e18 (KO)',
'WT vs KO (11)', 'WT vs KO (13)', 'WT vs KO (15)','WT vs KO (18)'
]
i = 0
values = []
for d in dfs:
sigD = d[d['padj'] < 0.05]
sigD = sigD[abs(sigD['log2FoldChange']) > logfc_cutoff]
num_sig = len(sigD)
values.append(num_sig)
print(f'{num_sig}, {labels[i]}')
i += 1
# Plot a simple bar chart.
barchart = Barchart(df, labels, values, title='Temporal genes anterior', order=labels)
barchart.palette=sns.color_palette([wt_colour, ko_colour, wt_colour, ko_colour, wt_colour, ko_colour,
wt_colour, ko_colour, wt_colour, ko_colour, wt_colour, ko_colour,
e11_colour, e13_colour, e15_colour, e18_colour])
barchart.plot()
save_fig(f'Barchart_anterior_time_numSigs')
u.dp(["Barchart used p.adj cutoff:", 0.05, "\nLogFC cutoff:", logfc_cutoff])
In [15]:
"""
---------------------------------------------------------------
Bar charts for time specific (anterior & posterior)
---------------------------------------------------------------
"""
dfs = [
df_dict[f'DEseq2_CNS_anterior_11_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_13_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_15_{date}.csv'],
df_dict[f'DEseq2_CNS_anterior_18_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_11_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_13_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_15_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_18_{date}.csv']
]
labels = [
'A. WT vs KO (11)', 'A. WT vs KO (13)', 'A. WT vs KO (15)','A. WT vs KO (18)',
'P. WT vs KO (11)', 'P. WT vs KO (13)', 'P. WT vs KO (15)','P. WT vs KO (18)'
]
i = 0
values = []
logfc_cutoff = 0.5
for d in dfs:
sigD = d[d['padj'] < 0.05]
sigD = sigD[abs(sigD['log2FoldChange']) > logfc_cutoff]
num_sig = len(sigD)
values.append(num_sig)
print(f'{num_sig}, {labels[i]}')
i += 1
# Plot a simple bar chart.
barchart = Barchart(dfs[0], labels, values, title='Temporal genes anterior', order=labels)
barchart.palette=sns.color_palette([e11_colour, e13_colour, e15_colour, e18_colour,
e11_colour, e13_colour, e15_colour, e18_colour])
barchart.plot()
save_fig('Barchart_anterior-posterior_time_numSigs')
u.dp(["Barchart used p.adj cutoff:", 0.05, "\nLogFC cutoff:", logfc_cutoff])
In [17]:
"""
---------------------------------------------------------------
Bar charts for time specific DE (posterior)
---------------------------------------------------------------
"""
dfs = [df_dict[f'DEseq2_CNS_posterior_wt_11-13_{date}.csv'], df_dict[f'DEseq2_CNS_posterior_ko_11-13_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_wt_13-15_{date}.csv'], df_dict[f'DEseq2_CNS_posterior_ko_13-15_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_wt_15-18_{date}.csv'], df_dict[f'DEseq2_CNS_posterior_ko_15-18_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_wt_11-15_{date}.csv'], df_dict[f'DEseq2_CNS_posterior_ko_11-15_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_wt_13-18_{date}.csv'], df_dict[f'DEseq2_CNS_posterior_ko_13-18_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_wt_11-18_{date}.csv'], df_dict[f'DEseq2_CNS_posterior_ko_11-18_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_11_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_13_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_15_{date}.csv'],
df_dict[f'DEseq2_CNS_posterior_18_{date}.csv']
]
labels = ['e11 vs e13 (WT)', 'e11 vs e13 (KO)',
'e13 vs e15 (WT)', 'e13 vs e15 (KO)',
'e15 vs e18 (WT)', 'e15 vs e18 (KO)',
'e11 vs e15 (WT)', 'e11 vs e15 (KO)',
'e13 vs e18 (WT)', 'e13 vs e18 (KO)',
'e11 vs e18 (WT)', 'e11 vs e18 (KO)',
'WT vs KO (11)', 'WT vs KO (13)', 'WT vs KO (15)','WT vs KO (18)'
]
i = 0
values = []
logfc_cutoff = 1.0
for d in dfs:
sigD = d[d['padj'] < 0.05]
sigD = sigD[abs(sigD['log2FoldChange']) > logfc_cutoff]
num_sig = len(sigD)
values.append(num_sig)
print(f'{num_sig}, {labels[i]}')
i += 1
# Plot a simple bar chart.
barchart = Barchart(df, labels, values, title='Temporal genes posterior', order=labels)
barchart.palette=sns.color_palette([wt_colour, ko_colour, wt_colour, ko_colour, wt_colour, ko_colour,
wt_colour, ko_colour, wt_colour, ko_colour, wt_colour, ko_colour,
e11_colour, e13_colour, e15_colour, e18_colour])
barchart.plot()
save_fig(f'Barchart_posterior_time_numSigs')
u.dp(["Barchart used p.adj cutoff:", 0.05, "\nLogFC cutoff:", logfc_cutoff])
8) Select log2(TMM + 1) as normalisation¶
We choose the normalisation to be the log2(TMM + 1) as we are keen to look across the genes rather than just run differential expression between the genes. As such, we need to load this add the information from the various differential expression analyses.
We also merge our WT vs KO experiments from before.
In [18]:
"""
--------------------------------------------------------
Choose normalisation method for data & add annoation
--------------------------------------------------------
"""
# Read in normalised RNAseq data
tmm = pd.read_csv(f'{r_dir}merged_df_FEATURE_COUNTS_tmm_{date}.csv')
tmm.rename(columns={ tmm.columns[0]: gene_id }, inplace = True)
for c in tmm:
if c != gene_id:
tmm[c] = np.log2(tmm[c].values + 1)
df_all = tmm.copy()
df_all[gene_id] = pd.to_numeric(df_all[gene_id])
# Read in each our our DE experiments from DEseq2
fb_df = df_dict[f'DEseq2_CNS_fb_{date}.csv']
mb_df = df_dict[f'DEseq2_CNS_mb_{date}.csv']
hb_df = df_dict[f'DEseq2_CNS_hb_{date}.csv']
sc_df = df_dict[f'DEseq2_CNS_sc_{date}.csv']
a11_df = df_dict[f'DEseq2_CNS_anterior_11_{date}.csv']
a13_df = df_dict[f'DEseq2_CNS_anterior_13_{date}.csv']
a15_df = df_dict[f'DEseq2_CNS_anterior_15_{date}.csv']
a18_df = df_dict[f'DEseq2_CNS_anterior_18_{date}.csv']
p11_df = df_dict[f'DEseq2_CNS_posterior_11_{date}.csv']
p13_df = df_dict[f'DEseq2_CNS_posterior_13_{date}.csv']
p15_df = df_dict[f'DEseq2_CNS_posterior_15_{date}.csv']
p18_df = df_dict[f'DEseq2_CNS_posterior_18_{date}.csv']
# Let's add all our significant results (and make sure we just indicate where they came from)
fb_df[gene_id] = pd.to_numeric(fb_df[gene_id])
df_all = df_all.merge(fb_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_fb'))
mb_df[gene_id] = pd.to_numeric(mb_df[gene_id])
df_all = df_all.merge(mb_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_mb'))
hb_df[gene_id] = pd.to_numeric(hb_df[gene_id])
df_all = df_all.merge(hb_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_hb'))
sc_df[gene_id] = pd.to_numeric(sc_df[gene_id])
df_all = df_all.merge(sc_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_sc'))
# Now we also want to add in the results from the anterior and posterior temporal changes
a11_df[gene_id] = pd.to_numeric(a11_df[gene_id])
df_all = df_all.merge(a11_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_a11'))
a13_df[gene_id] = pd.to_numeric(a13_df[gene_id])
df_all = df_all.merge(a13_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_a13'))
a15_df[gene_id] = pd.to_numeric(a15_df[gene_id])
df_all = df_all.merge(a15_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_a15'))
a18_df[gene_id] = pd.to_numeric(a18_df[gene_id])
df_all = df_all.merge(a18_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_a18'))
p11_df[gene_id] = pd.to_numeric(p11_df[gene_id])
df_all = df_all.merge(p11_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_p11'))
p13_df[gene_id] = pd.to_numeric(p13_df[gene_id])
df_all = df_all.merge(p13_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_p13'))
p15_df[gene_id] = pd.to_numeric(p15_df[gene_id])
df_all = df_all.merge(p15_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_p15'))
p18_df[gene_id] = pd.to_numeric(p18_df[gene_id])
df_all = df_all.merge(p18_df, how='left', left_on=gene_id,
right_on=gene_id, suffixes=('', '_p18'))
# Rename the padj from forebrain and logfc
df_all.rename(columns={'padj': 'padj_fb', 'log2FoldChange': 'log2FoldChange_fb'}, inplace=True)
# Drop columns that don't have gene names --> leave this out for now.
#df_all = df_all.dropna(subset=[gene_name])
# This will be on p values that weren't in there
p_cols = [c for c in df_all.columns if 'padj' in c]
# Make them NS i.e. given them a p value of 1.0
df_all[p_cols] = df_all[p_cols].fillna(value=1.0)
# Fill everything else with 0's so we don't have NaNs in our counts etc
df_all = df_all.fillna(value=0)
df_all = df_all.reset_index() # Ensure the index is reset each time we rejoin everything
u.dp(['Length of dataframe', len(df_all)])
# We don't want duplicate gene IDs so we drop columns with duplicates in these values
no_na_merged = df_all.drop_duplicates()
# Collect entrez IDs
ens_ids = annot_df['ensembl_gene_id'].values
chr_ids = annot_df['chromosome_name'].values
gene_id_to_ens = {}
gene_id_to_chr = {}
for i, g in enumerate(annot_df[gene_id].values):
if not gene_id_to_ens.get(g):
gene_id_to_ens[g] = ens_ids[i]
if not gene_id_to_chr.get(g):
gene_id_to_chr[g] = chr_ids[i]
ens_ids = []
chrs = []
for g in no_na_merged[gene_id].values:
ens_ids.append(gene_id_to_ens.get(g))
chrs.append(gene_id_to_chr.get(g))
# Add ensembl gene id
no_na_merged['ensembl_gene_id'] = ens_ids
no_na_merged['chrs'] = chrs
# Organise columns a bit nicer
cols = [gene_id, gene_name, 'ensembl_gene_id', 'chrs',
'wt11fb1', 'wt11fb2', 'wt13fb1', 'wt13fb2', 'wt15fb1', 'wt15fb2', 'wt18fb1', 'wt18fb2',
'wt11mb1', 'wt11mb2', 'wt13mb1', 'wt13mb2', 'wt15mb1', 'wt15mb2', 'wt18mb1', 'wt18mb2',
'wt11hb1', 'wt11hb2', 'wt13hb1', 'wt13hb2', 'wt15hb1', 'wt15hb2', 'wt18hb1', 'wt18hb2',
'wt11sc1', 'wt11sc2', 'wt13sc1', 'wt13sc2', 'wt15sc1', 'wt15sc2', 'wt18sc1', 'wt18sc2',
'ko11fb1', 'ko11fb2', 'ko13fb1', 'ko13fb2', 'ko15fb1', 'ko15fb2', 'ko18fb1', 'ko18fb2',
'ko11mb1', 'ko11mb2', 'ko13mb1', 'ko13mb2', 'ko15mb1', 'ko15mb2', 'ko18mb1', 'ko18mb2',
'ko11hb1', 'ko11hb2', 'ko13hb1', 'ko13hb2', 'ko15hb1', 'ko15hb2', 'ko18hb1', 'ko18hb2',
'ko11sc1', 'ko11sc2', 'ko13sc1', 'ko13sc2', 'ko15sc1', 'ko15sc2', 'ko18sc1', 'ko18sc2',
]
for c in no_na_merged:
if c not in cols:
cols.append(c)
df_all = no_na_merged[cols]
df_all = df_all.reset_index()
# Ensure the gene name hasn't been messed up
gene_names = []
pseudo_id = []
for g in df_all[gene_id].values:
gene_n = gene_id_to_name.get(g)
gene_names.append(gene_n)
pseudo_id.append(f'{g}-{gene_n}')
df_all[gene_name] = gene_names
u.dp(['Length of dataframe (no dups reset index)', len(df_all)])
Testing gender of samples¶
Here we just confirm the gender of the samples based on X and Y reads.
In [22]:
y_chr = df_all[df_all['chrs'] == 'Y']
x_chr = df_all[df_all['chrs'] == 'X']
wt_cols = [c for c in df_all.columns if 'merged' not in c and 'wt' in c]
for w in wt_cols:
plt.hist(x_chr[w].values)
plt.title(f'WT X: {w}')
plt.show()
wt_cols = [c for c in df_all.columns if 'merged' not in c and 'wt' in c]
for w in wt_cols:
plt.hist(y_chr[w].values)
plt.title(f'WT Y: {w}')
plt.show()
ko_cols = [c for c in df_all.columns if 'merged' not in c and 'ko' in c]
for w in ko_cols:
plt.hist(x_chr[w].values)
plt.title(f'KO X: {w}')
plt.show()
for w in ko_cols:
plt.hist(y_chr[w].values)
plt.title(f'KO Y: {w}')
plt.show()
9) Visualise the merged data top most significant genes¶
Here we visualise the top genes by padj in the forebrain and validate that these align with the other results.
We observe the Hox genes as the most diverged.
In [23]:
from matplotlib.colors import ListedColormap
from sklearn.decomposition import PCA
from sklearn.preprocessing import MinMaxScaler
from sciviso import Scatterplot
import umap
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from sciviso import Vis
# Plot PCA of the forebrain genes i.e. so we see how time separates
def get_cond_colour(c):
if 'ko' in c:
return ko_colour
elif 'wt' in c:
return wt_colour
return '#FFFFFF'
def plt_heatmap(df, idxs, title='', plt_sig=True, mark=None, row_colours=None, cond=''):
# Now let's look at a couple of these, and how similar these are to the other data
heatmap_df = pd.DataFrame()
# Select our genes of interest
if plt_sig:
ps = df['padj'].values[idxs]
p_idxs = np.where(ps < 0.05)
heatmap_df[gene_name] = df[gene_name].values[idxs][p_idxs]
else:
heatmap_df[gene_name] = df[gene_name].values[idxs]
cols = []
col_colours = []
cond_colours = []
tissue_colours = []
time_colours = []
for c in df.columns:
if mark == None:
if (cond in c) and 'merged' in c:
cname = c.split('_')[0]
print(c)
if plt_sig:
heatmap_df[cname] = df[c].values[idxs][p_idxs]
else:
heatmap_df[cname] = df[c].values[idxs]
cols.append(cname)
time_colours.append(get_time_colour(c))
tissue_colours.append(get_tissue_colour(c))
cond_colours.append(get_cond_colour(c))
else:
if mark in c and hist_metric in c:
if plt_sig:
heatmap_df[c] = np.log2(df[c].values[idxs][p_idxs] + 1)
else:
heatmap_df[c] = np.log2(df[c].values[idxs] + 1)
time_colours.append(get_time_colour(c))
tissue_colours.append(get_tissue_colour(c))
cond_colours.append(get_mark_colour(c))
cols.append(c)
col_colours = [cond_colours, tissue_colours, time_colours]
heatmap = Heatmap(heatmap_df, cols, gene_name, linewidths=0.0, cmap="Purples", x_tick_labels=0,
figsize=(2, 1.5), vmin=0, vmax=9,
title=f'{title}', cluster_cols=False, cluster_rows=True,
col_colours=col_colours, row_colours=row_colours)
ax = heatmap.plot()
ax.ax_heatmap.tick_params(direction='out', length=2, width=0.5)
ax.ax_heatmap.tick_params(labelsize=6)
ax.ax_heatmap.tick_params(axis='x', which='major', pad=2.0)
ax.ax_heatmap.tick_params(axis='y', which='major', pad=2.0)
pplot()
save_fig(title)
plt.show()
Plot PCA with each of the features¶
In [33]:
cols, edges = [], []
e13, e13_colours = [], []
e18, e18_colours = [], []
e15, e15_colours = [], []
tissue_colours = []
time_colours = []
cond_colours = []
wt_tissue_colours = []
wt_time_colours = []
wt_cond_colours = []
ko_tissue_colours = []
ko_time_colours = []
ko_cond_colours = []
c_i = 0
wt_pts = []
ko_pts = []
line_styles = []
for c in df_all.columns:
if ('wt' in c or 'ko' in c) and '11' not in c:
cols.append(c)
if 'wt' in c:
wt_tissue_colours.append(get_tissue_colour(c))
wt_time_colours.append(get_time_colour(c))
wt_pts.append(c_i)
elif 'ko' in c:
ko_tissue_colours.append(get_tissue_colour(c))
ko_time_colours.append(get_time_colour(c))
ko_pts.append(c_i)
if '13' in c:
e13.append(c_i)
e13_colours.append(get_tissue_colour(c))
tissue_colours.append(get_tissue_colour(c))
time_colours.append(get_time_colour(c))
cond_colours.append(get_cond_colour(c))
elif '15' in c:
e15.append(c_i)
e15_colours.append(get_tissue_colour(c))
tissue_colours.append(get_tissue_colour(c))
time_colours.append(get_time_colour(c))
cond_colours.append(get_cond_colour(c))
elif '18' in c:
e18.append(c_i)
e18_colours.append(get_tissue_colour(c))
# Make the edges the colour of our condition
edges.append(get_cond_colour(c))
tissue_colours.append(get_tissue_colour(c))
time_colours.append(get_time_colour(c))
cond_colours.append(get_cond_colour(c))
if 'wt' in c:
line_styles.append('-')
else:
line_styles.append('-')
c_i += 1
vals = (df_all[cols].values).T #np.log2(df_all[cols].values + 1).T
fb_pca = PCA(n_components=2)
fb_pca_values = fb_pca.fit_transform(vals)
var_ratio = fb_pca.fit(vals).explained_variance_ratio_
"""
---------------------------------------------------------------
Plot PCA
---------------------------------------------------------------
"""
plt.rcParams['figure.figsize'] = [2, 2]
plt.scatter(fb_pca_values[e13,0], fb_pca_values[e13,1], linestyle=line_styles, c=e13_colours, s=100, marker=">", edgecolors=edges, linewidths=1.5)
plt.scatter(fb_pca_values[e15,0], fb_pca_values[e15,1], linestyle=line_styles, c=e15_colours, s=100, marker="o", edgecolors=edges, linewidths=1.5)
plt.scatter(fb_pca_values[e18,0], fb_pca_values[e18,1], linestyle=line_styles, c=e18_colours, s=100, marker="X", edgecolors=edges, linewidths=1.5)
plt.title(f'PCA VAR: 0: {var_ratio[0]}, 1: {var_ratio[1]}')
save_fig(f'PCA_ne11_fb')
# Now we want to fit everything except the gene IDs
plt.show()
"""
---------------------------------------------------------------
PCA coloured by tissue
---------------------------------------------------------------
"""
plt.rcParams['figure.figsize'] = [2, 2]
plt.scatter(fb_pca_values[wt_pts,0], fb_pca_values[wt_pts,1], linestyle=line_styles, c=wt_tissue_colours, s=100, marker="o", edgecolors='black', linewidths=0.5)
plt.scatter(fb_pca_values[ko_pts,0], fb_pca_values[ko_pts,1], linestyle=line_styles, c=ko_tissue_colours, s=100, marker="X", edgecolors='black', linewidths=0.5)
plt.title(f'PCA VAR: 0: {var_ratio[0]}, 1: {var_ratio[1]}')
save_fig(f'PCA_ne11_tissue')
# Now we want to fit everything except the gene IDs
plt.show()
"""
---------------------------------------------------------------
PCA coloured by time
---------------------------------------------------------------
"""
plt.rcParams['figure.figsize'] = [2, 2]
plt.scatter(fb_pca_values[wt_pts,0], fb_pca_values[wt_pts,1], linestyle=line_styles, c=wt_time_colours, s=100, marker="o", edgecolors='black', linewidths=0.5)
plt.scatter(fb_pca_values[ko_pts,0], fb_pca_values[ko_pts,1], linestyle=line_styles, c=ko_time_colours, s=100, marker="X", edgecolors='black', linewidths=0.5)
plt.title(f'PCA VAR: 0: {var_ratio[0]}, 1: {var_ratio[1]}')
save_fig(f'PCA_ne11_time')
# Now we want to fit everything except the gene IDs
plt.show()
Merge samples and do sample clustermap¶
Here we just do a sample clustermap
In [88]:
"""
---------------------------------------------------------------
Sample clustermap for anterior late stage
---------------------------------------------------------------
"""
# Smooth out the columns in the data frame i.e. for the clones we only put in the mean of the two replicates
log2_df = pd.DataFrame()
cols_to_merge = [c for c in df_all.columns if 'wt' in c or 'ko' in c]
col_names = []
col_values = []
i = 0
tissue_colours = []
cond_colours = []
time_colours = []
while(i < len(cols_to_merge)):
if ('mb' in cols_to_merge[i] or 'fb' in cols_to_merge[i]) and ('15' in cols_to_merge[i] or '18' in cols_to_merge[i]):
log2_df[cols_to_merge[i][:-1]] = 0.5 * ((df_all[cols_to_merge[i]].values + 1) +
df_all[cols_to_merge[i + 1]].values)
tissue_colours.append(get_tissue_colour(cols_to_merge[i]))
time_colours.append(get_time_colour(cols_to_merge[i]))
cond_colours.append(get_cond_colour(cols_to_merge[i]))
print("merged", cols_to_merge[i], cols_to_merge[i + 1])
i += 2
plt.rcParams['figure.figsize'] = [3, 3]
row_colors_t = [cond_colours, tissue_colours, time_colours]
corr = log2_df.corr()
sns.clustermap(corr,
xticklabels=corr.columns.values,
yticklabels=corr.columns.values, cmap='RdBu_r', row_cluster=True,
col_cluster=True, row_colors=row_colors_t)
save_fig(f'Heatmap_sample_cluster')
plt.show()
"""
---------------------------------------------------------------
Sample clustermap for all
---------------------------------------------------------------
"""
# Smooth out the columns in the data frame i.e. for the clones we only put in the mean of the two replicates
log2_df = pd.DataFrame()
cols_to_merge = [c for c in df_all.columns if 'wt' in c or 'ko' in c]
col_names = []
col_values = []
i = 0
tissue_colours = []
cond_colours = []
time_colours = []
while(i < len(cols_to_merge)):
log2_df[cols_to_merge[i][:-1]] = 0.5 * ((df_all[cols_to_merge[i]].values + 1) +
df_all[cols_to_merge[i + 1]].values)
tissue_colours.append(get_tissue_colour(cols_to_merge[i]))
time_colours.append(get_time_colour(cols_to_merge[i]))
cond_colours.append(get_cond_colour(cols_to_merge[i]))
print("merged", cols_to_merge[i], cols_to_merge[i + 1])
i += 2
plt.rcParams['figure.figsize'] = [3, 3]
row_colors_t = [cond_colours, tissue_colours, time_colours]
corr = log2_df.corr()
sns.clustermap(corr,
xticklabels=corr.columns.values,
yticklabels=corr.columns.values, cmap='RdBu_r', row_cluster=True,
col_cluster=True, row_colors=row_colors_t)
save_fig(f'Heatmap_sample_cluster_all')
"""
---------------------------------------------------------------
Sample clustermap for E11
---------------------------------------------------------------
"""
# Smooth out the columns in the data frame i.e. for the clones we only put in the mean of the two replicates
log2_df = pd.DataFrame()
cols_to_merge = [c for c in df_all.columns if 'wt' in c or 'ko' in c]
col_names = []
col_values = []
i = 0
tissue_colours = []
cond_colours = []
time_colours = []
while(i < len(cols_to_merge)):
if '11' in cols_to_merge[i]:
log2_df[cols_to_merge[i][:-1]] = 0.5 * ((df_all[cols_to_merge[i]].values + 1) +
df_all[cols_to_merge[i + 1]].values)
tissue_colours.append(get_tissue_colour(cols_to_merge[i]))
time_colours.append(get_time_colour(cols_to_merge[i]))
cond_colours.append(get_cond_colour(cols_to_merge[i]))
print("merged", cols_to_merge[i], cols_to_merge[i + 1])
i += 2
plt.rcParams['figure.figsize'] = [3, 3]
row_colors_t = [cond_colours, tissue_colours, time_colours]
corr = log2_df.corr()
sns.clustermap(corr,
xticklabels=corr.columns.values,
yticklabels=corr.columns.values, cmap='RdBu_r', row_cluster=True,
col_cluster=True, row_colors=row_colors_t)
save_fig(f'Heatmap_sample_cluster_Ell')
"""
---------------------------------------------------------------
Sample clustermap for E13
---------------------------------------------------------------
"""
# Smooth out the columns in the data frame i.e. for the clones we only put in the mean of the two replicates
log2_df = pd.DataFrame()
cols_to_merge = [c for c in df_all.columns if 'wt' in c or 'ko' in c]
col_names = []
col_values = []
i = 0
tissue_colours = []
cond_colours = []
time_colours = []
while(i < len(cols_to_merge)):
if '13' in cols_to_merge[i]:
log2_df[cols_to_merge[i][:-1]] = 0.5 * ((df_all[cols_to_merge[i]].values + 1) +
df_all[cols_to_merge[i + 1]].values)
tissue_colours.append(get_tissue_colour(cols_to_merge[i]))
time_colours.append(get_time_colour(cols_to_merge[i]))
cond_colours.append(get_cond_colour(cols_to_merge[i]))
print("merged", cols_to_merge[i], cols_to_merge[i + 1])
i += 2
plt.rcParams['figure.figsize'] = [3, 3]
row_colors_t = [cond_colours, tissue_colours, time_colours]
corr = log2_df.corr()
sns.clustermap(corr,
xticklabels=corr.columns.values,
yticklabels=corr.columns.values, cmap='RdBu_r', row_cluster=True,
col_cluster=True, row_colors=row_colors_t)
save_fig(f'Heatmap_sample_cluster_El3')
"""
---------------------------------------------------------------
Sample clustermap for E15
---------------------------------------------------------------
"""
# Smooth out the columns in the data frame i.e. for the clones we only put in the mean of the two replicates
log2_df = pd.DataFrame()
cols_to_merge = [c for c in df_all.columns if 'wt' in c or 'ko' in c]
col_names = []
col_values = []
i = 0
tissue_colours = []
cond_colours = []
time_colours = []
while(i < len(cols_to_merge)):
if '15' in cols_to_merge[i]:
log2_df[cols_to_merge[i][:-1]] = 0.5 * ((df_all[cols_to_merge[i]].values + 1) +
df_all[cols_to_merge[i + 1]].values)
tissue_colours.append(get_tissue_colour(cols_to_merge[i]))
time_colours.append(get_time_colour(cols_to_merge[i]))
cond_colours.append(get_cond_colour(cols_to_merge[i]))
print("merged", cols_to_merge[i], cols_to_merge[i + 1])
i += 2
plt.rcParams['figure.figsize'] = [3, 3]
row_colors_t = [cond_colours, tissue_colours, time_colours]
corr = log2_df.corr()
sns.clustermap(corr,
xticklabels=corr.columns.values,
yticklabels=corr.columns.values, cmap='RdBu_r', row_cluster=True,
col_cluster=True, row_colors=row_colors_t)
save_fig(f'Heatmap_sample_cluster_El5')
"""
---------------------------------------------------------------
Sample clustermap for E18
---------------------------------------------------------------
"""
# Smooth out the columns in the data frame i.e. for the clones we only put in the mean of the two replicates
log2_df = pd.DataFrame()
cols_to_merge = [c for c in df_all.columns if 'wt' in c or 'ko' in c]
col_names = []
col_values = []
i = 0
tissue_colours = []
cond_colours = []
time_colours = []
while(i < len(cols_to_merge)):
if '15' in cols_to_merge[i]:
log2_df[cols_to_merge[i][:-1]] = 0.5 * ((df_all[cols_to_merge[i]].values + 1) +
df_all[cols_to_merge[i + 1]].values)
tissue_colours.append(get_tissue_colour(cols_to_merge[i]))
time_colours.append(get_time_colour(cols_to_merge[i]))
cond_colours.append(get_cond_colour(cols_to_merge[i]))
print("merged", cols_to_merge[i], cols_to_merge[i + 1])
i += 2
plt.rcParams['figure.figsize'] = [3, 3]
row_colors_t = [cond_colours, tissue_colours, time_colours]
corr = log2_df.corr()
sns.clustermap(corr,
xticklabels=corr.columns.values,
yticklabels=corr.columns.values, cmap='RdBu_r', row_cluster=True,
col_cluster=True, row_colors=row_colors_t)
save_fig(f'Heatmap_sample_cluster_El8')
Plot heatmaps of samples RNAseq¶
Plot the top and bottom of the heatmaps of each of the DE analyses.
In [34]:
"""
---------------------------------------------------------------
Select top 20 DE genes (by padj in Forebrain)
---------------------------------------------------------------
"""
n_genes = 10
n_g_neg = -10
# Select the locfc that is the forebrain WT vs KO
merged_df = df_all.copy()
cols_to_merge = [c for c in merged_df.columns if 'wt' in c or 'ko' in c]
col_names = []
col_values = []
values = merged_df[cols_to_merge].values
i = 0
while(i < len(cols_to_merge)):
merged_df[f'{cols_to_merge[i][:-1]}_merged-rep'] = 0.5 * (merged_df[cols_to_merge[i]].values +
merged_df[cols_to_merge[i + 1]].values)
print("merged", cols_to_merge[i], cols_to_merge[i + 1])
i += 2
log2FoldChange = merged_df['stat'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_genes)[:n_genes]) # np.argpartition(x, -k)[-k:]
col_colours = ['#6666ff'] * n_genes + ['#990033'] * n_genes
"""
---------------------------------------------------------------
Plot heatmap
---------------------------------------------------------------
"""
plt_heatmap(merged_df, selected_idxs, f'Top genes', False, row_colours=None, cond="")
plt.show()
selected_idxs = list(np.argpartition(log2FoldChange, n_g_neg)[n_g_neg:])
plt_heatmap(merged_df, selected_idxs, f'Bottom genes', False, row_colours=None, cond="")
plt.show()
"""
---------------------------------------------------------------
Setup PCA
---------------------------------------------------------------
"""
cols, edges = [], []
e13, e13_colours = [], []
e18, e18_colours = [], []
e15, e15_colours = [], []
c_i = 0
line_styles = []
for c in df_all.columns:
if ('wt' in c or 'ko' in c) and '11' not in c:
cols.append(c)
if '13' in c:
e13.append(c_i)
e13_colours.append(get_tissue_colour(c))
elif '15' in c:
e15.append(c_i)
e15_colours.append(get_tissue_colour(c))
elif '18' in c:
e18.append(c_i)
e18_colours.append(get_tissue_colour(c))
# Make the edges the colour of our condition
edges.append(get_cond_colour(c))
if 'wt' in c:
line_styles.append('-')
else:
line_styles.append('-')
c_i += 1
vals = np.log2(df_all[cols].values + 1).T
fb_pca = PCA(n_components=2)
fb_pca_values = fb_pca.fit_transform(vals)
var_ratio = fb_pca.fit(vals).explained_variance_ratio_
"""
---------------------------------------------------------------
Plot PCA
---------------------------------------------------------------
"""
plt.rcParams['figure.figsize'] = [2, 2]
plt.scatter(fb_pca_values[e13,0], fb_pca_values[e13,1], linestyle=line_styles, c=e13_colours, s=100, marker=">", edgecolors=edges, linewidths=1.5)
plt.scatter(fb_pca_values[e15,0], fb_pca_values[e15,1], linestyle=line_styles, c=e15_colours, s=100, marker="o", edgecolors=edges, linewidths=1.5)
plt.scatter(fb_pca_values[e18,0], fb_pca_values[e18,1], linestyle=line_styles, c=e18_colours, s=100, marker="X", edgecolors=edges, linewidths=1.5)
plt.title(f'PCA VAR: 0: {var_ratio[0]}, 1: {var_ratio[1]}')
save_fig(f'PCA_ne11_fb')
# Now we want to fit everything except the gene IDs
plt.show()
In [35]:
"""
---------------------------------------------------------------
Select top 20 DE genes (by padj in Forebrain)
---------------------------------------------------------------
"""
n_genes = 10
n_g_neg = -10
val = 'stat'
# Select the locfc that is the forebrain WT vs KO
for cond in ['a11', 'a13', 'a15', 'a18', 'p11', 'p13', 'p15', 'p18']:
e18_gene_df = merged_df.copy()
log2FoldChange = e18_gene_df[f'{val}_{cond}'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_genes)[:n_genes]) # np.argpartition(x, -k)[-k:]
col_colours = ['#6666ff'] * n_genes + ['#990033'] * n_genes
"""
---------------------------------------------------------------
Plot heatmap
---------------------------------------------------------------
"""
plt_heatmap(merged_df, selected_idxs, f'Downregulated genes {cond}', False, row_colours=None)
plt.show()
log2FoldChange = e18_gene_df[f'log2FoldChange_{cond}'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_g_neg)[n_g_neg:])
plt_heatmap(merged_df, selected_idxs, f'Upregulated genes {cond}', False, row_colours=None)
plt.show()
In [36]:
"""
---------------------------------------------------------------
Select top 20 DE genes (by padj in Forebrain)
---------------------------------------------------------------
"""
n_genes = 10
n_g_neg = -10
# Select the locfc that is the forebrain WT vs KO
e18_gene_df = merged_df.copy()
log2FoldChange = e18_gene_df['stat_p18'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_genes)[:n_genes]) # np.argpartition(x, -k)[-k:]
col_colours = ['#6666ff'] * n_genes + ['#990033'] * n_genes
"""
---------------------------------------------------------------
Plot heatmap
---------------------------------------------------------------
"""
plt_heatmap(merged_df, selected_idxs, f'Downregulated genes HB-SC e18', False, row_colours=None)
plt.show()
log2FoldChange = e18_gene_df['log2FoldChange_p18'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_g_neg)[n_g_neg:])
plt_heatmap(merged_df, selected_idxs, f'Upregulated genes HB-SC e18', False, row_colours=None)
plt.show()
In [37]:
"""
---------------------------------------------------------------
Select top 20 DE genes (by padj in Forebrain)
---------------------------------------------------------------
"""
n_genes = 10
n_g_neg = -10
# Select the locfc that is the forebrain WT vs KO
e18_gene_df = merged_df.copy()
log2FoldChange = e18_gene_df['log2FoldChange_sc'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_genes)[:n_genes]) # np.argpartition(x, -k)[-k:]
col_colours = ['#6666ff'] * n_genes + ['#990033'] * n_genes
"""
---------------------------------------------------------------
Plot heatmap
---------------------------------------------------------------
"""
plt_heatmap(merged_df, selected_idxs, f'Downregulated genes SC', False, row_colours=None)
plt.show()
log2FoldChange = e18_gene_df['log2FoldChange_sc'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_g_neg)[n_g_neg:])
plt_heatmap(merged_df, selected_idxs, f'Upregulaed genes SC', False, row_colours=None)
plt.show()
In [38]:
"""
---------------------------------------------------------------
Select top 20 DE genes (by padj in Forebrain)
---------------------------------------------------------------
"""
n_genes = 10
n_g_neg = -10
# Select the locfc that is the forebrain WT vs KO
e18_gene_df = merged_df.copy()
log2FoldChange = e18_gene_df['log2FoldChange_hb'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_genes)[:n_genes]) # np.argpartition(x, -k)[-k:]
col_colours = ['#6666ff'] * n_genes + ['#990033'] * n_genes
"""
---------------------------------------------------------------
Plot heatmap
---------------------------------------------------------------
"""
plt_heatmap(merged_df, selected_idxs, f'Downregulated genes HB', False, row_colours=None)
plt.show()
log2FoldChange = e18_gene_df['log2FoldChange_hb'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_g_neg)[n_g_neg:])
plt_heatmap(merged_df, selected_idxs, f'Upregulaed genes HB', False, row_colours=None)
plt.show()
In [39]:
"""
---------------------------------------------------------------
Select top 20 DE genes (by padj in Forebrain)
---------------------------------------------------------------
"""
n_genes = 10
n_g_neg = -10
# Select the locfc that is the forebrain WT vs KO
e18_gene_df = merged_df.copy()
log2FoldChange = e18_gene_df['log2FoldChange_mb'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_genes)[:n_genes]) # np.argpartition(x, -k)[-k:]
col_colours = ['#6666ff'] * n_genes + ['#990033'] * n_genes
"""
---------------------------------------------------------------
Plot heatmap
---------------------------------------------------------------
"""
plt_heatmap(merged_df, selected_idxs, f'Downregulated genes MB', False, row_colours=None)
plt.show()
log2FoldChange = e18_gene_df['log2FoldChange_mb'].values
selected_idxs = list(np.argpartition(log2FoldChange, n_g_neg)[n_g_neg:])
plt_heatmap(merged_df, selected_idxs, f'Upregulaed genes MB', False, row_colours=None)
plt.show()
Heatmaps of marker genes in WT and KO¶
In [40]:
gene_markers_sep = [['Foxg1', 'Tbr1', 'Emx1', 'Eomes'],
['Bhlhe23', 'En2', 'Fev', 'Phox2b'],
['Hoxb9', 'Hoxc9', 'Hoxd9', 'Hoxc12'],
['Hoxd8', 'Hoxd9', 'Hoxd10', 'Hoxd11', 'Hoxd12', 'Hoxd13', 'Hoxa7', 'Hoxa9', 'Hoxa10', 'Hoxa11',
'Hoxa13', 'Hoxb9', 'Hoxb13', 'Hoxc8', 'Hoxc9', 'Hoxc10', 'Hoxc11', 'Hoxc12', 'Hoxc13'],
['Ccna1', 'Ccna2', 'Ccnd1', 'Ccnd2', 'Ccnd3', 'Ccne1', 'Ccne2', 'Cdc25a',
'Cdc25b', 'Cdc25c', 'E2f1', 'E2f2', 'E2f3', 'Mcm10', 'Mcm5', 'Mcm3', 'Mcm2', 'Cip2a'],
['Cdkn1a', 'Cdkn1b', 'Cdkn1c', 'Cdkn2a', 'Cdkn2b', 'Cdkn2c', 'Cdkn2d'],
['Sox1', 'Sox2', 'Sox3'],
['Snap25', 'Syt1', 'Slc32a1','Slc17a6', 'Syn1'],
['Cspg4', 'Aqp4', 'Slc6a11', 'Olig1', 'Igfbp3']
]
marker_labels_sep = ['Forebrain', 'Midbrain', 'Hindbrain', 'Spinalcord',
'Proliferation', 'Neg. reg. Cell Cycle',
'Progenitors', 'Neurons', 'Glia'
]
In [41]:
for i, genes in enumerate(gene_markers_sep):
selected_idxs = []
for j, gene in enumerate(df_all[gene_name].values):
if gene in genes:
selected_idxs.append(j)
plt_heatmap(merged_df, selected_idxs, marker_labels_sep[i] + " genes WT", False, row_colours=None, cond='wt')
plt_heatmap(merged_df, selected_idxs, marker_labels_sep[i] + " genes KO", False, row_colours=None, cond='ko')
plt.show()
10) Add Epigenetic data from Encode¶
We downloaded and processed the data according to the notebook DownloadEncode. Now we want to merge the peak files with our RNAseq data.
Here we keep the signal with the highest value as multiple may have been added.
Keeping in mind, we have potentially the same peak assigned to many genes. This is particularly evident in the H3K27me3 broad peak data.
First we download an annotation from We use sciloc2gene v1.0.0 to annotate this information.
Note you may need to check the reference used for the peak files (this was not needed since we also used Ensembl to map our gene's and find their location).
Need to download the annotation from: https://www.encodeproject.org/references/ENCSR425FOI/ ENCFF871VGROpen file informationDownload gtf genome reference mm10 M21
Convert the GTF to a bed and then sort the bed file (convert using gtf2bed v0.11.0 downloaded frmo https://github.com/fls-bioinformatics-core/GFFUtils/tree/master) gtf2bed ENCFF871VGR.gtf > ENCFF871VGR_mm10-m21_encode.bed
Then we convert to a CSV file and sort it by TSS.
def run_parallel(self):
data_dir = 'data/encode/histone_modifications/beds_cns/encode/sorted_bed/'
files = os.listdir(data_dir)
files_to_run = []
for filename in files:
files_to_run.append(filename)
# Run in paralell
pool = ThreadPool(12)
results = pool.map(run_bed, files_to_run)
def run_bed(filename):
base_dir = 'data/encode/histone_modifications/beds_cns/encode/'
data_dir = f'{base_dir}sorted_bed/'
output_dir = f'{base_dir}scie2g_28102020/'
overlap_method = 'in_promoter'
if 'H3K9me3' in filename or 'H3K36me3' in filename:
overlap_method = 'overlaps'
#ensembl_gene_id, entrezgene_id, external_gene_name, chromosome_name, start_position, end_position, strand
bed = Bed(f'{data_dir}{filename}', overlap_method=overlap_method,
output_bed_file=f'{output_dir}selected_peaks/{filename}',
buffer_after_tss=1500,
buffer_before_tss=5000,
buffer_gene_overlap=1500,
chr_idx=0, start_idx=1,
end_idx=2, peak_value=6, header_extra='8,9',
header="chr,start,end,signal,qvalue,peak",
gene_start=4, gene_end=5, gene_direction=6, gene_chr=3, gene_name=0
)
# Use your annotation file i.e. for humans
bed.set_annotation_from_file('../supps/mm10Sorted_mmusculus_gene_ensembl-GRCm38.p6.csv')
# Now we can run the assign values
bed.assign_locations_to_genes()
bed.save_loc_to_csv(f'{output_dir}parsed_files/{filename[:-3]}csv', keep_unassigned=True)In [42]:
"""
--------------------------------------------------------
Add Epigenetic data to the RNAseq dataframe
--------------------------------------------------------
"""
%matplotlib inline
def add_peak_file(df, filename):
# Make a dictionary from the key being the gene id to the value (width)
mapping = {}
mapping_signals = {}
file_df = pd.read_csv(f'{chip_dir}{filename}')
# Keep the max value
signals = file_df['signal'].values
widths = file_df['width'].values
qvalues = file_df['qvalue'].values # q-value at position 9 i.e. 9 - 1 = 8
# Only keep peaks that were significant i.e. qvalue > -1 * log10(0.05) ~ 1.3
for i, g in enumerate(file_df[gene_id].values):
qvalue = qvalues[i]
if qvalue > 1.3:
if not mapping.get(g):
mapping[g] = widths[i]
mapping_signals[g] = signals[i]
else:
mapping[g] = widths[i] if widths[i] > mapping[g] else mapping[g]
# Keep the signal that corresponds to the widest peak
mapping_signals[g] = signals[i] if widths[i] > mapping[g] else mapping_signals[g]
# Now we want to iterate through the genes in the df and add them
widths_ordered = []
signals_ordered = []
fname = filename[:-4]
for g in df[gene_id].values:
widths_ordered.append(mapping.get(g))
signals_ordered.append(mapping_signals.get(g))
df[f'{fname}_width'] = widths_ordered
df[f'{fname}_signal'] = signals_ordered
return df
# Now we want to add our chipseq data from Encode
chip_dir = os.path.join(input_dir, f'scie2g_b2500_g500/parsed_files/')
files = os.listdir(chip_dir)
chip_files = []
mark_locations = []
for f in files:
if '.csv' in f and '.swp' not in f and 'embryonic_' in f:
chip_files.append(f)
chip_files.sort()
# Add them to the dataframe
dfs = []
for filename in chip_files:
u.dp(["Adding", filename])
df_all = add_peak_file(df_all, filename)
11) Filter the merged dataframe¶
Since we have merged so many dataframes there are likely to be many missing values.
We want to first drop any rows that are duplicates & 2) fill in null values with 0's (otherwise we'll get errors when running analyses).
We also want to organise the columns of the dataframe (this makes it easier for charts later on).
Lastly, we are only interested in genes that had a significant change between WT and KO (we want to desipher WHY this occured), so we drop any data with a non-sig result.
Filtering steps:
1) Keep genes if padj < 0.05 in at least one experiment, THEN
2) Keep genes if genes were expressed in WT or KO tissues with at least an average TMM of 2, THEN
3) Keep genes if genes at least one experiment recorded a log fold change of > 1.0
This dataset will be used for our VAE analysis.
In [43]:
# Keep only the significant data i.e. the data where at least one of the p values was < 0.05
genes_to_keep = []
ps = ((1.0 * (df_all['padj_sc'].values < 0.05)) + (1.0 * (df_all['padj_mb'].values < 0.05))
+ (1.0 * (df_all['padj_hb'].values < 0.05)) + (1.0 * (df_all['padj_fb'].values < 0.05))
+ (1.0 * (df_all['padj_a11'].values < 0.05))
+ (1.0 * (df_all['padj_a13'].values < 0.05))
+ (1.0 * (df_all['padj_a15'].values < 0.05))
+ (1.0 * (df_all['padj_a18'].values < 0.05)))
df_all['p_mask'] = ps
# Now we also want to make sure that in each of these datasets, we have at least 2 TMM average in either WT or KO
# We don't want the genes that have super low expression since that is not interesting for our downstream analyses
wt_cols = [c for c in df_all.columns if 'wt' in c]
ko_cols = [c for c in df_all.columns if 'ko' in c]
tmm_cutoff = 0.5
expression_mask = ((1.0 * (np.mean(df_all[wt_cols].values, axis=1) >= tmm_cutoff)) + (1.0 * (np.mean(df_all[ko_cols].values, axis=1) >= tmm_cutoff)))
df_all['expression_mask'] = expression_mask
# Make a DF that only has genes that are significant at the 0.05 level
df_sig = df_all[df_all['p_mask'] > 0]
# Make a DF that only has genes that are minimally expressed in at least one condition
df_sig = df_sig[df_sig['expression_mask'] > 0]
u.dp(['Number of significant rows: ', len(df_sig)])
# Let's also make sure that we can keep have genes that have at least an absolute logfc of 0.5
log_str = 'log2FoldChange'
logfc_cols = [c for c in df_sig.columns if log_str in c]
mean_logfc = np.mean(abs(df_sig[logfc_cols].values), axis=1)
# df_sig['mean_logfc'] = mean_logfc
# df_sig = df_sig[df_sig['mean_logfc'] > 0.5]
cutoff = 1.0
logfc = ((1.0 * (abs(df_sig[f'{log_str}_sc'].values) > cutoff))
+ (1.0 * (abs(df_sig[f'{log_str}_mb'].values) > cutoff))
+ (1.0 * (abs(df_sig[f'{log_str}_hb'].values) > cutoff))
+ (1.0 * (abs(df_sig[f'{log_str}_fb'].values) > cutoff))
+ (1.0 * (abs(df_sig[f'{log_str}_a11'].values) > cutoff))
+ (1.0 * (abs(df_sig[f'{log_str}_a13'].values) > cutoff))
+ (1.0 * (abs(df_sig[f'{log_str}_a15'].values) > cutoff))
+ (1.0 * (abs(df_sig[f'{log_str}_a18'].values) > cutoff))
)
df_sig['logfc_mask'] = logfc
df_sig_2 = df_sig[df_sig['logfc_mask'] > 2] # Require two of the conditions it to be met
df_sig_1 = df_sig[df_sig['logfc_mask'] > 0] # Require two of the conditions it to be met
u.dp(['Number of significant rows with an absolute FC > 1.0: ', len(df_sig_1)])
u.dp(['Number of significant rows with 2 values absolute FC > 2.0: ', len(df_sig_2)])
# Save all to dataframes
df_sig_1.to_csv(f'{output_dir}df-sig_1_epi-2500_{date}.csv', index=False)
df_sig_2.to_csv(f'{output_dir}df-consistent_epi-2500_{date}.csv', index=False)
df_sig.to_csv(f'{output_dir}df-significant_epi-2500_{date}.csv', index=False)
df_all.to_csv(f'{output_dir}df-all_epi-2500_{date}.csv', index=False)
# --------------------------------------------------------------------------------
# Number of significant rows: 12532
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------
# Number of significant rows with an absolute FC > 1.0: 2944
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------
# Number of significant rows with 2 values absolute FC > 2.0: 1928
# --------------------------------------------------------------------------------
gene_markers_sep = [['Emx1', 'Eomes', 'Tbr1', 'Foxg1', 'Lhx6'],
['En1', 'En2', 'Lmx1a', 'Bhlhe23', 'Sall4'],
['Hoxb1', 'Krox20', 'Fev', 'Hoxd3', 'Phox2b'],
['Hoxd8', 'Hoxd9', 'Hoxd10', 'Hoxd11', 'Hoxd12', 'Hoxd13', 'Hoxa7', 'Hoxa9', 'Hoxa10', 'Hoxa11',
'Hoxa13', 'Hoxb9', 'Hoxb13', 'Hoxc8', 'Hoxc9', 'Hoxc10', 'Hoxc11', 'Hoxc12', 'Hoxc13'],
['Ccna1', 'Ccna2', 'Ccnd1', 'Ccnd2', 'Ccnd3', 'Ccne1', 'Ccne2', 'Cdc25a',
'Cdc25b', 'Cdc25c', 'E2f1', 'E2f2', 'E2f3', 'Mcm10', 'Mcm5', 'Mcm3', 'Mcm2', 'Cip2a'],
['Cdkn1a', 'Cdkn1b', 'Cdkn1c', 'Cdkn2a', 'Cdkn2b', 'Cdkn2c', 'Cdkn2d'],
['Sox2', 'Sox1', 'Sox3', 'Hes1', 'Hes5'],
['Snap25', 'Syt1', 'Slc32a1','Slc17a6', 'Syn1'],
['Cspg4', 'Aqp4', 'Slc6a11', 'Olig1', 'Igfbp3'],
['Foxg1'],
['En2'],
['Phox2b'],
['Hoxc9'],
['Sox3']
]
for ml in gene_markers_sep:
print("\n")
for m in ml:
for g in df_sig_1[gene_name].values:
try:
if g == m:
print(g)
except:
x = 0
12) Setup functions for simple analyses¶
We want to test whether some simple things are enriched or not. i.e. do we see different annotations if we look at the genes that were 1) significant a) with H3K27me3, b) without H3K27me3, 2) non-significant a) with H3K27me3, b) without H3K27me3.
For each of these groups we want to:
1) test for annotation enrichment
2) plot the overal HM profile
3) plot the overal gene expression profile.In [44]:
import venn
from matplotlib.colors import ListedColormap
def plot_venn(gene_sets, labels, title, save=True, show=True, colours=None):
output = f'venn_{title.replace(" ", "-")}'
dset = {}
for i, l in enumerate(labels):
dset[l] = gene_sets[i]
if len(gene_sets) > 3:
cmap = ListedColormap(sns.color_palette(colours))
venn.venn(dset, cmap=cmap, fontsize=6, figsize=(2, 2))
else:
venn3(gene_sets, set_labels=labels)
if save:
save_fig(output)
if show:
plt.show()
In [45]:
from statsmodels.stats.multitest import multipletests
import string
from scipy import stats
from sciviso import Heatmap
import matplotlib
def run_annot_plot(df_bg, df_fg, title=""):
changes = []
order = ['Pr-A', 'Pr-W', 'Pr-B', 'Pr-F', 'En-Sd', 'En-Sp',
'En-W', 'En-Pd', 'En-Pp', 'Tr-S', 'Tr-P', 'Tr-I',
'Hc-P', 'Hc-H', 'NS']
total = len(df_bg) * 1.0
total_sig = len(df_fg) * 1.0
perc_all = []
perc_sig = []
for o in order:
o_all = len(df_bg[df_bg['peak_value'] == o])
o_sig = len(df_fg[df_fg['peak_value'] == o])
perc_all.append(o_all)
perc_sig.append(o_sig)
if o_all == 0 and o_sig == 0:
changes.append(0)
else:
changes.append(((o_sig) - (o_all)))
odds_ratios = []
pvalues = []
for i, o in enumerate(order):
# Do a FET on each one
oddsratio, pvalue = stats.fisher_exact([[perc_sig[i], perc_all[i]],
[total_sig - perc_sig[i],
total - perc_all[i]]])
print(o, oddsratio, pvalue,[[perc_sig[i], perc_all[i]],
[total_sig - perc_sig[i],
total - perc_all[i]]])
odds_ratios.append(oddsratio)
pvalues.append(pvalue)
reg, padj, a, b = multipletests(pvalues, alpha=0.1,
method='fdr_bh', returnsorted=False)
p_sigs = []
for p in padj:
if p > 0.05:
p_sigs.append('')
elif p <= 0.05 and p > 0.01:
p_sigs.append('*')
elif p <= 0.01 and p > 0.001:
p_sigs.append('**')
elif p <= 0.001 and p > 0.0001:
p_sigs.append('***')
elif p <= 0.0001:
p_sigs.append('****')
else:
print(p)
p_sigs.append('')
fig, ax = plt.subplots(figsize=(1.5,1.0))
c = [grey] * 14
c[2] = "green"
c[12] = "green"
plt.bar(order, odds_ratios, color=c, linewidth=0.5, edgecolor='black')
rects = ax.patches
# Make some labels.
labels = [f'{p_sigs[i]}' for i in range(0, len(rects))]
for rect, label in zip(rects, labels):
height = rect.get_height()
ax.text(rect.get_x() + rect.get_width() / 2, height, label,
ha='center', va='bottom')
plt.ylim(0, 7)
plt.yticks(np.arange(0, 7, 2.0))
ax.set_title(f'Odds ratio FET for significant dataset {title} compared to all genes ChromHMM annotations')
ax.set_xlabel('', fontsize=6)
ax.set_ylabel('Odds ratio', fontsize=6)
ax.set_xticklabels(order, rotation=90, ha="center")
ax.tick_params(direction='out', length=2, width=0.5)
ax.spines['bottom'].set_linewidth(0.5)
ax.spines['top'].set_linewidth(0)
ax.spines['left'].set_linewidth(0.5)
ax.spines['right'].set_linewidth(0)
ax.tick_params(labelsize=6)
ax.tick_params(axis='x', which='major', pad=0)
ax.tick_params(axis='y', which='major', pad=0)
save_fig(f'OddsRatio-ChromHMM-{title}')
plt.show()
tissues = ['forebrain', 'midbrain', 'hindbrain', 'neural-tube',
'embryonic-facial-prominence', 'limb',
'heart', 'liver']
marks = ['H3K9me3', 'H3K36me3', 'H3K27me3', 'H3K27ac', 'H3K4me1', 'H3K4me2', 'H3K4me3','H3K9ac']
def get_median_naninc(df, mark, hist_metric="", time_1="", time_2="", time_3="", tissue="brain"):
cols = []
for c in df.columns:
if tissue in c and mark in c and hist_metric in c and (time_1 in c or time_2 in c or time_3 in c):
cols.append(c)
# get nan median
vals = np.nanmean(df[cols].values, axis=1)
return vals
def plot_mark_heatmap(df, idxs, title, mark_cutoff=3.0):
mark_df = pd.DataFrame()
mark_values = []
mean_col = []
num_genes = 0
if idxs is None:
num_genes = len(df)
else:
num_genes = len(idxs) * 1.0
for t in [['10.5', '11.5', '12.5'], ['15.5', '15.5', '16.5']]:
tissue_titles = []
marks_title = []
for m in marks:
mark_col = []
for tissue in tissues:
median_all_data = get_median_naninc(df, m, "signal", t[0], t[1], t[2], tissue)
if idxs is None:
median_genes = median_all_data
else:
median_genes = median_all_data[idxs]
has_mark = 1.0 * len(np.where(median_genes > mark_cutoff)[0])
if has_mark == 0:
mark_col.append(0)
else:
mark_col.append(has_mark/num_genes) #np.nan_to_num(np.nanmean(median_genes)))
if string.capwords(tissue) not in tissue_titles:
if tissue == 'embryonic-facial-prominence':
if 'E.F.P' not in tissue_titles:
tissue_titles.append('E.F.P')
else:
tissue_titles.append(string.capwords(tissue))
marks_title.append(string.capwords(m))
mark_df[string.capwords(m)] = mark_col
mark_df['Tissue'] = tissue_titles
heatmap = Heatmap(mark_df, marks_title, 'Tissue', vmin=0, vmax=1, cmap='Greens', figsize=(2, 2),
title=f'{title} {"-".join(t)}', cluster_rows=False, cluster_cols=False)
heatmap.plot()
pplot()
print(mark_df.head())
save_fig(f'mark_all_tissues_signal-{title}_{"-".join(t)}')
plt.show()
13) Setup the groups¶
1) Significnat (p < 0.05 in at least one experiment (anterior time or spatial wt/ko))
1.a) has H3K27me3 mark (signal at e16.5 > 3.0)
1.a.i) Expressed (mean (log2(TMM + 1) > 2 in WT OR KO)
1.a.j) Not expressed
1.b) doesn't have H3K27me3 (signal at e16.5 < 3.0)
1.b.i) Expressed
1.b.j) Not expressed
2) Not significant
2.a) Has (as above)
2.a.i) Exprerssed
2.a.j) Not expressed
2.b) Has not (as above)
2.b.i) Expressed
2.b.j) Not expressedIn [68]:
tmm_cutoff = 0.5
df_bg = df_all.copy() # Make a copy so we don't override the other one
# Also have a look at how the marks plot over genes that are not significant
chrom_dir = os.path.join(input_dir, "supps", "annot")
chromhmm_annot = pd.read_csv(f'{chrom_dir}/e16.5_forebrain_15_segments.csv')
print(len(chromhmm_annot), len(df_bg))
# Make sure we only keep 1 annotation per gene
chromhmm_annot = chromhmm_annot.groupby('external_gene_name').first()
# Merge this with our dataframe
df_bg = df_bg.merge(chromhmm_annot, on='external_gene_name', how='left', suffixes=('', '_chmm'))
# Ensure the new length is the same as the old length
ps = ((1.0 * (df_all['padj_sc'].values < 0.05)) + (1.0 * (df_all['padj_mb'].values < 0.05))
+ (1.0 * (df_all['padj_hb'].values < 0.05)) + (1.0 * (df_all['padj_fb'].values < 0.05))
+ (1.0 * (df_all['padj_a11'].values < 0.05))
+ (1.0 * (df_all['padj_a13'].values < 0.05))
+ (1.0 * (df_all['padj_a15'].values < 0.05))
+ (1.0 * (df_all['padj_a18'].values < 0.05)))
df_bg['p_mask'] = ps
df_bg = df_bg.fillna(0)
ko_cols = [c for c in df_bg if 'ko' in c]
wt_cols = [c for c in df_bg if 'wt' in c]
expression_mask = ((1.0 * (np.mean(df_all[wt_cols].values, axis=1) >= tmm_cutoff)) + (1.0 * (np.mean(df_all[ko_cols].values, axis=1) >= tmm_cutoff)))
df_bg['expression_mask'] = expression_mask
df_bg = df_bg[df_bg['expression_mask'] > tmm_cutoff]
cutoff = 1.0
log_str = 'log2FoldChange'
logfc = ((1.0 * (abs(df_bg[f'{log_str}_sc'].values) > cutoff))
+ (1.0 * (abs(df_bg[f'{log_str}_mb'].values) > cutoff))
+ (1.0 * (abs(df_bg[f'{log_str}_hb'].values) > cutoff))
+ (1.0 * (abs(df_bg[f'{log_str}_fb'].values) > cutoff))
+ (1.0 * (abs(df_bg[f'{log_str}_a11'].values) > cutoff))
+ (1.0 * (abs(df_bg[f'{log_str}_a13'].values) > cutoff))
+ (1.0 * (abs(df_bg[f'{log_str}_a15'].values) > cutoff))
+ (1.0 * (abs(df_bg[f'{log_str}_a18'].values) > cutoff))
)
df_bg['logfc_mask'] = logfc
# Add in the H3K mark --> we'll just look in the brain regions at e16.5
h3k_cols = [c for c in df_bg.columns if 'brain' in c and '16.5' in c and 'signal' in c and 'H3K27me3' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
ectopic = ((1.0 * (np.mean(df_bg[wt_cols].values, axis=1) <= tmm_cutoff)) + (1.0 * (np.mean(df_bg[ko_cols].values, axis=1) >= tmm_cutoff)))
df_bg['ectopic'] = ectopic
df_sig = df_bg[df_bg['p_mask'] > 0]
df_ns = df_bg[df_bg['p_mask'] == 0]
h3k_cutoff = 0.0
df_sig_hasH3 = df_sig[df_sig['hk3_mean'] > h3k_cutoff]
df_hasH3_exp = df_bg[df_bg['hk3_mean'] > h3k_cutoff]
df_noH3_exp = df_bg[df_bg['hk3_mean'] == 0]
df_noH3_exp = df_noH3_exp[df_noH3_exp['expression_mask'] > tmm_cutoff]
pert_cutoff = 2
df_sig_hasH3_pert = df_sig_hasH3[df_sig_hasH3['logfc_mask'] > pert_cutoff] # Consistently perturbed
df_sig_hasH3_unpert = df_sig_hasH3[df_sig_hasH3['logfc_mask'] <= pert_cutoff] # Inconsistently perturbed
df_sig_hasH3_ectopic = df_sig_hasH3[df_sig_hasH3['ectopic'] == pert_cutoff]
df_sig_hasH3_exp = df_sig_hasH3[df_sig_hasH3['expression_mask'] > tmm_cutoff]
df_sig_hasH3_exp_pert = df_sig_hasH3_exp[df_sig_hasH3_exp['logfc_mask'] > pert_cutoff] # Consistently perturbed
df_sig_hasH3_exp_unpert = df_sig_hasH3_exp[df_sig_hasH3_exp['logfc_mask'] <= pert_cutoff] # Inonsistently perturbed
df_sig_hasH3_not_exp = df_sig_hasH3[df_sig_hasH3['expression_mask'] <= tmm_cutoff]
df_sig_noH3 = df_sig[df_sig['hk3_mean'] <= h3k_cutoff]
df_sig_noH3_pert = df_sig_noH3[df_sig_noH3['logfc_mask'] > pert_cutoff] # Consistently perturbed
df_sig_noH3_unpert = df_sig_noH3[df_sig_noH3['logfc_mask'] <= pert_cutoff] # Inonsistently perturbed
df_sig_noH3_exp = df_sig_noH3[df_sig_noH3['expression_mask'] > tmm_cutoff]
# #df_bg = df_bg[df_bg['logfc_mask'] > 2]
df_sig_noH3_exp_pert = df_sig_noH3_exp[df_sig_noH3_exp['logfc_mask'] > pert_cutoff]
df_sig_noH3_exp_unpert = df_sig_noH3_exp[df_sig_noH3_exp['logfc_mask'] <= pert_cutoff]
df_sig_noH3_ectopic = df_sig_noH3[df_sig_noH3['ectopic'] == 2]
df_sig_noH3_notexp = df_sig_noH3[df_sig_noH3['expression_mask'] <= tmm_cutoff]
df_ns_hasH3 = df_ns[df_ns['hk3_mean'] > h3k_cutoff]
df_ns_hasH3_exp = df_ns_hasH3[df_ns_hasH3['expression_mask'] > tmm_cutoff]
df_ns_hasH3_not_exp = df_ns_hasH3[df_ns_hasH3['expression_mask'] <= tmm_cutoff]
df_ns_noH3 = df_ns[df_ns['hk3_mean'] <= h3k_cutoff]
df_ns_noH3_exp = df_ns_noH3[df_ns_noH3['expression_mask'] > tmm_cutoff]
df_ns_noH3_notexp = df_ns_noH3[df_ns_noH3['expression_mask'] <= tmm_cutoff]
df_ns_noH3_ectopic = df_ns_noH3[df_ns_noH3['ectopic'] == 2]
df_ns_H3_ectopic = df_ns_hasH3[df_ns_hasH3['ectopic'] == 2]
dfs = [df_bg, df_sig, df_ns, df_sig_hasH3, df_ns_hasH3, df_sig_noH3, df_ns_noH3,
df_sig_hasH3_pert, df_sig_hasH3_unpert, df_sig_noH3_pert, df_sig_noH3_unpert,
df_sig_hasH3_exp_pert, df_sig_hasH3_exp_unpert, df_sig_noH3_exp_pert, df_sig_noH3_exp_unpert,
df_sig_hasH3_exp, df_sig_hasH3_not_exp, df_sig_noH3_exp, df_sig_noH3_notexp,
df_ns_hasH3_exp, df_ns_hasH3_not_exp, df_ns_noH3_exp, df_ns_noH3_notexp,
df_sig_hasH3_ectopic, df_sig_noH3_ectopic]
df_dict = {'All genes': df_bg,
'Genes with H3K27me3': df_hasH3_exp,
'Genes without H3K27me3': df_noH3_exp,
'Genes with at least 1 sig': df_sig,
'Non-sig genes': df_ns,
'Sig with H3K27me3': df_sig_hasH3,
'NS with H3K27me3': df_ns_hasH3,
'Sig unmarked': df_sig_noH3,
'NS unmarked': df_ns_noH3,
'Sig Marked Perturbed': df_sig_hasH3_pert,
'Sig Marked un-Perturbed': df_sig_hasH3_unpert,
'Sig Unmarked Perturbed': df_sig_noH3_pert,
'Sig Unmarked un-Perturbed': df_sig_noH3_unpert,
'Exp. Sig Marked Perturbed': df_sig_hasH3_exp_pert,
'Exp. Sig Marked un-Perturbed': df_sig_hasH3_exp_unpert,
'Exp. Sig Unmarked Perturbed': df_sig_noH3_exp_pert,
'Exp. Sig Unmarked un-Perturbed': df_sig_noH3_exp_unpert,
'Sig genes with H3K27me3 and expression': df_sig_hasH3_exp,
'Sig genes with H3K27me3 NO expression': df_sig_hasH3_not_exp,
'Sig genes unmarked and expression': df_sig_noH3_exp,
'Sig genes unmarked NO expression': df_sig_noH3_notexp,
'NS genes with H3K27me3 and expression': df_ns_hasH3_exp,
'NS genes with H3K27me3 NO expression': df_ns_hasH3_not_exp,
'NS genes unmarked and expression': df_ns_noH3_exp,
'NS genes unmarked NO expression': df_ns_noH3_notexp,
'Sig ectopic marked': df_sig_hasH3_ectopic,
'Sig ectopic unmarked': df_sig_noH3_ectopic,
'NS ectopic marked': df_ns_H3_ectopic,
'NS ectopic unmarked': df_ns_noH3_ectopic,
}
emg = [g for g in df_dict['Sig ectopic marked'][gene_name].values]
print(", ".join(emg))
emg = [g for g in df_dict['Sig ectopic unmarked'][gene_name].values]
print(", ".join(emg))
Plot correlation between FB logFC and median H3K27me3 mark¶
In [51]:
%matplotlib inline
from scipy.stats import spearmanr
from sciviso import Histogram, Scatterplot
from numpy.polynomial.polynomial import polyfit
import statsmodels.api as sm
# Do this on the consistently affected dataset
df_const = pd.read_csv(f'{output_dir}df-consistent_epi-2500_{date}.csv')
#calculate Spearman Rank correlation and corre
# Also do the main markers
fb_genes = ['Emx1', 'Eomes', 'Tbr1', 'Foxg1', 'Lhx6']
mb_genes = ['En1', 'En2', 'Lmx1a', 'Bhlhe23', 'Sall4']
hb_genes = ['Phox2b', 'Krox20', 'Fev', 'Hoxb1', 'Hoxd3']
sc_genes = ['Hoxd8', 'Hoxd9', 'Hoxd10', 'Hoxd11','Hoxa7', 'Hoxa9', 'Hoxa10',
'Hoxb9', 'Hoxb13', 'Hoxc8', 'Hoxc9', 'Hoxc10', 'Hoxc11', 'Hoxc12', 'Hoxc13']
# Fill NAs in genes that didn't have the H3K27me3 mark
df_const_na = df_const.fillna(0)
h3k_cols = [c for c in df_const_na.columns if 'brain' in c and 'H3K27me3' in c and 'signal' in c]
Y = np.nan_to_num(np.nanmedian(df_const_na[h3k_cols].values, axis=1))
df_const_na['Median_H3K'] = Y
X = df_const_na['log2FoldChange_fb'].values
results = sm.OLS(Y,sm.add_constant(X)).fit()
print(results.summary())
corr, _ = spearmanr(X, Y)
u.dp(["H3K27me3 median signal in brain regions vs FB logFC:\n", 'Spearmans correlation: %.3f' % corr])
# now run the scatter
s = Scatterplot(df_const_na,'Median_H3K', 'log2FoldChange_fb', title='H3K27me3 vs FB logFC', colour="green",
points_to_annotate=['Foxg1', 'Otx2', 'Pax2', 'Hoxd9'],
annotation_label=gene_name, add_legend=False)
s.opacity=0.1
ax = s.plot()
ax.set_xlim(-3, 20)
ax.set_ylim(-5, 12)
save_fig(f'Scatter_H3K27me3vsLogFC')
plt.show()
"""
--------------------------------------------------------
Do the same with H3K27ac as a control
--------------------------------------------------------
"""
h3k_cols = [c for c in df_const_na.columns if 'brain' in c and 'H3K27ac' in c and 'signal' in c]
Y = np.nan_to_num(np.nanmedian(df_const_na[h3k_cols].values, axis=1))
df_const_na['Median_H3K'] = Y
X = df_const_na['log2FoldChange_fb'].values
results = sm.OLS(Y,sm.add_constant(X)).fit()
print(results.summary())
corr, _ = spearmanr(X, Y)
u.dp(["H3K27ac median signal in brain regions vs FB logFC:\n", 'Spearmans correlation: %.3f' % corr])
# now run the scatter
s = Scatterplot(df_const_na,'Median_H3K', 'log2FoldChange_fb', title='H3K27ac vs FB logFC', colour="grey",
points_to_annotate=['Foxg1', 'Otx2', 'Pax2', 'Hoxd9'],
annotation_label=gene_name, add_legend=False)
s.opacity=0.1
ax = s.plot()
ax.set_xlim(-3, 20)
ax.set_ylim(-5, 12)
save_fig(f'Scatter_H3K27acvsLogFC')
In [52]:
h3k_cols = [c for c in df_bg.columns if 'brain' in c and '10.5' in c and 'signal' in c and 'H3K27me3' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
e10 = df_bg[df_bg['hk3_mean'] > h3k_cutoff][gene_name].values
h3k_cols = [c for c in df_bg.columns if 'brain' in c and '11.5' in c and 'signal' in c and 'H3K27me3' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
e11 = df_bg[df_bg['hk3_mean'] > h3k_cutoff][gene_name].values
h3k_cols = [c for c in df_bg.columns if 'brain' in c and '12.5' in c and 'signal' in c and 'H3K27me3' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
e12 = df_bg[df_bg['hk3_mean'] > h3k_cutoff][gene_name].values
h3k_cols = [c for c in df_bg.columns if 'brain' in c and '13.5' in c and 'signal' in c and 'H3K27me3' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
e13 = df_bg[df_bg['hk3_mean'] > h3k_cutoff][gene_name].values
h3k_cols = [c for c in df_bg.columns if 'brain' in c and '14.5' in c and 'signal' in c and 'H3K27me3' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
e14 = df_bg[df_bg['hk3_mean'] > h3k_cutoff][gene_name].values
h3k_cols = [c for c in df_bg.columns if 'brain' in c and '15.5' in c and 'signal' in c and 'H3K27me3' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
e15 = df_bg[df_bg['hk3_mean'] > h3k_cutoff][gene_name].values
h3k_cols = [c for c in df_bg.columns if 'brain' in c and '16.5' in c and 'signal' in c and 'H3K27me3' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
e16 = df_bg[df_bg['hk3_mean'] > h3k_cutoff][gene_name].values
### -----------
plt.plot([10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5],
[len(e10), len(e11), len(e12), len(e13), len(e14), len(e15), len(e16)]
)
afct = set(df_sig[gene_name].values)
cafct = set(df_sig[df_sig['logfc_mask'] > pert_cutoff][gene_name].values)
plt.plot([10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5],
[len(set(afct) & set(e10)), len(set(afct) & set(e11)),
len(set(afct) & set(e12)), len(set(afct) & set(e13)), len(set(afct) & set(e14)),
len(set(afct) & set(e15)), len(set(afct) & set(e16))]
)
plt.plot([10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5],
[len(set(cafct) & set(e10)), len(set(cafct) & set(e11)),
len(set(cafct) & set(e12)), len(set(cafct) & set(e13)), len(set(cafct) & set(e14)),
len(set(cafct) & set(e15)), len(set(cafct) & set(e16))]
)
plt.show()
# ==============================
plt.plot([10.5, 12.5, 13.5, 14.5, 15.5, 16.5],
[len(e10), len(e12), len(e13), len(e14), len(e15), len(e16)], label="marked"
)
plt.plot([10.5, 12.5, 13.5, 14.5, 15.5, 16.5],
[len(df_bg) - len(e10), len(df_bg) - len(e12), len(df_bg) - len(e13), len(df_bg) - len(e14), len(df_bg) - len(e15), len(df_bg) - len(e16)]
, label="bg - marked")
afct = set(df_sig[gene_name].values)
cafct = set(df_sig[df_sig['logfc_mask'] > pert_cutoff][gene_name].values)
afct = set(df_sig[gene_name].values)
plt.plot([10.5, 12.5, 13.5, 14.5, 15.5, 16.5],
[len(set(afct) & set(e10)),
len(set(afct) & set(e12)), len(set(afct) & set(e13)), len(set(afct) & set(e14)),
len(set(afct) & set(e15)), len(set(afct) & set(e16))], label="Intersect, affect and marked"
)
plt.plot([10.5, 12.5, 13.5, 14.5, 15.5, 16.5],
[len(set(afct)) - len(set(afct) & set(e10)),
len(set(afct)) - len(set(afct) & set(e12)), len(set(afct)) - len(set(afct) & set(e13)),
len(set(afct)) - len(set(afct) & set(e14)),
len(set(afct)) - len(set(afct) & set(e15)), len(set(afct)) - len(set(afct) & set(e16))],
label="unmarked affected"
)
plt.plot([10.5, 12.5, 13.5, 14.5, 15.5, 16.5],
[len(set(cafct) & set(e10)),
len(set(cafct) & set(e12)), len(set(cafct) & set(e13)), len(set(cafct) & set(e14)),
len(set(cafct) & set(e15)), len(set(cafct) & set(e16))],
label="consist. affected, marked"
)
plt.legend(loc="upper left")
plt.show()
#------------
h3k_cols = [c for c in df_bg.columns if 'brain' in c and 'signal' in c and 'H3K27me3' in c and '16.5' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
plot_venn((set(df_bg[df_bg['hk3_mean'] == 0][gene_name].values),
set(df_sig[gene_name].values), set(df_sig[df_sig['logfc_mask'] > pert_cutoff][gene_name].values)),
('Un-Marked genes', 'Affected Genes', 'Const. Affected Genes'), 'venn_unmaked_marked_const_aff')
#------------
h3k_cols = [c for c in df_bg.columns if 'brain' in c and 'signal' in c and 'H3K27me3' in c and '10.5' in c]
df_bg['hk3_mean'] = np.nanmean(df_bg[h3k_cols].values, axis=1)
plot_venn((set(df_bg[df_bg['hk3_mean'] > h3k_cutoff][gene_name].values),
set(df_sig[gene_name].values), set(df_sig[df_sig['logfc_mask'] > pert_cutoff][gene_name].values)),
('Marked genes', 'Affected Genes', 'Const. Affected Genes'), 'venn_marked_const_aff')
In [53]:
# Print out the size of each DF
for label, d in df_dict.items():
u.dp([label, len(d)])
13) Test for enrichment of specific marks using the Gorkin et al annotations¶
In [69]:
# Save each DF so we can test each of them for ORA (just save the gene name to avoid dups.)
for label, d in df_dict.items():
print(label)
dn = d[gene_id]
dn.to_csv(os.path.join(output_dir, f'{label.replace(" ", "_")}_gene-name_{date}.csv'), index=False)
In [56]:
"""
--------------------------------------------------------
Test for enriched gorkin et al states
--------------------------------------------------------
"""
for label, d in df_dict.items():
run_annot_plot(df_bg, d, title=f'Annot {label}')
In [58]:
for label, d in df_dict.items():
try:
logFC_up = d[d['log2FoldChange_fb'] > 0.5]
logFC_down = d[d['log2FoldChange_fb'] < -0.5]
u.dp([label, len(d), 100 * (len(d)/len(df_all))])
print('LogFC up: ', len(logFC_up), 100 * (len(logFC_up)/len(d)))
print('LogFC down: ', len(logFC_down), 100 * (len(logFC_down)/len(d)))
except:
print("unable to run:", label)
15) Plot the general trend of the RNA seq¶
In [59]:
fb_genes = ['Emx1', 'Eomes', 'Tbr1', 'Foxg1', 'Lhx6']
mb_genes = ['En1', 'En2', 'Lmx1a', 'Bhlhe23', 'Sall4']
hb_genes = ['Phox2b', 'Krox20', 'Fev', 'Hoxb1', 'Hoxd3']
sc_genes = ['Hoxd8', 'Hoxd9', 'Hoxd10', 'Hoxd11','Hoxa7', 'Hoxa9', 'Hoxa10',
'Hoxb9', 'Hoxb13', 'Hoxc8', 'Hoxc9', 'Hoxc10', 'Hoxc11', 'Hoxc12', 'Hoxc13']
progenitors = ['Sox2', 'Sox1', 'Sox3', 'Hes1', 'Hes5']
neurons = ['Tubb3', 'Snap25', 'Syt1', 'Slc32a1','Slc17a6']
glia = ['Pdgfra', 'Cspg4', 'Aqp4', 'Egfr', 'Slc6a11']
from collections import defaultdict
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from statannot import add_stat_annotation
hox_genes = sc_genes
forebrain_genes = fb_genes
h3cols = []
fb_rna_cols = []
sc_rna_cols = []
# Do it in this way so they are correctly ordered
for c in df_bg.columns:
if 'H3K27me' in c and 'signal' in c:
h3cols.append(c)
elif 'wt' in c:
if 'fb' in c:
fb_rna_cols.append(c)
if 'sc' in c:
sc_rna_cols.append(c)
for c in df_bg.columns:
if 'ko' in c:
if 'fb' in c:
fb_rna_cols.append(c)
if 'sc' in c:
sc_rna_cols.append(c)
wt_colour = "#AADFF1"
def plot_gene_boxplot(df, title, cluster_id, cols):
boxplot = Boxplot(df, gene_name, cols[0])
boxplot.label_font_size = 6
boxplot.title_font_size = 8
sns.set(rc={'figure.figsize': (1.5, 1.5), 'font.family': 'sans-serif',
'font.sans-serif': 'Arial', 'font.size': 8.0}, style='ticks')
boxplot.palette = sci_colour
box_df = boxplot.format_data_for_boxplot(df, cols, gene_name, df[gene_name].values)
is_hox = []
sns.set_style("ticks")
for c in box_df['Samples'].values:
if 'Hox' in c:
is_hox.append('Hox')
else:
is_hox.append('Hox like')
box_df['is_hox'] = is_hox
boxplot = Boxplot(box_df, "Conditions", "Values", add_stats=False, add_dots=False,
order=cols)
boxplot.label_font_size = 6
boxplot.title_font_size = 8
plt.rcParams['figure.figsize'] = (1.5,1.5)
sns.set(rc={'figure.figsize': (1.5, 1.5), 'font.family': 'sans-serif',
'font.sans-serif': 'Arial', 'font.size': 8.0}, style='ticks')
boxplot.palette = sci_colour
ax = boxplot.plot()
ax.tick_params(labelsize=6)
plt.title(title, {'fontsize': 8, 'fontweight': 700})
ax.set_ylim(0, 2.5)
c = 0
for b in ax.artists:
if c < len(cols)/2:
b.set_facecolor(wt_colour)
else:
b.set_facecolor(ko_colour)
c += 1
save_fig(f'{title}')
plt.show()
avgs_fb = [['wt11fb1',
'wt11fb2'],
['wt13fb1',
'wt13fb2'],
['wt15fb1',
'wt15fb2'],
['wt18fb1',
'wt18fb2'],
['ko11fb1',
'ko11fb2'],
['ko13fb1',
'ko13fb2'],
['ko15fb1',
'ko15fb2'],
['ko18fb1',
'ko18fb2']]
avgs_sc = [['wt11sc1',
'wt11sc2'],
['wt13sc1',
'wt13sc2'],
['wt15sc1',
'wt15sc2'],
['wt18sc1',
'wt18sc2'],
['ko11sc1',
'ko11sc2'],
['ko13sc1',
'ko13sc2'],
['ko15sc1',
'ko15sc2'],
['ko18sc1',
'ko18sc2']]
avgs_hb = [['wt11hb1',
'wt11hb2'],
['wt13hb1',
'wt13hb2'],
['wt15hb1',
'wt15hb2'],
['wt18hb1',
'wt18hb2'],
['ko11hb1',
'ko11hb2'],
['ko13hb1',
'ko13hb2'],
['ko15hb1',
'ko15hb2'],
['ko18hb1',
'ko18hb2']]
avgs_mb = [['wt11mb1',
'wt11mb2'],
['wt13mb1',
'wt13mb2'],
['wt15mb1',
'wt15mb2'],
['wt18mb1',
'wt18mb2'],
['ko11mb1',
'ko11mb2'],
['ko13mb1',
'ko13mb2'],
['ko15mb1',
'ko15mb2'],
['ko18mb1',
'ko18mb2']]
for label, df in df_dict.items():
if 'ectopic' in label:
avgs_df = pd.DataFrame()
avgs_df[gene_name] = df[gene_name].values
fb_cols = []
for f in avgs_fb:
new_col = f'{f[0][:2]} {f[0][2:4]} FB'
avgs_df[new_col] = (df[f[0]].values + df[f[1]].values) / 2.0
fb_cols.append(new_col)
sc_cols = []
for f in avgs_sc:
new_col = f'{f[0][:2]} {f[0][2:4]} SC'
avgs_df[new_col] = (df[f[0]].values + df[f[1]].values) / 2.0
sc_cols.append(new_col)
mb_cols = []
for f in avgs_mb:
new_col = f'{f[0][:2]} {f[0][2:4]} MB'
avgs_df[new_col] = (df[f[0]].values + df[f[1]].values) / 2.0
mb_cols.append(new_col)
hb_cols = []
for f in avgs_hb:
new_col = f'{f[0][:2]} {f[0][2:4]} HB'
avgs_df[new_col] = (df[f[0]].values + df[f[1]].values) / 2.0
hb_cols.append(new_col)
plot_gene_boxplot(avgs_df, f'{label} FB', '', fb_cols)
plot_gene_boxplot(avgs_df, f'{label} MB', '', mb_cols)
plot_gene_boxplot(avgs_df, f'{label} HB', '', hb_cols)
plot_gene_boxplot(avgs_df, f'{label} SC', '', sc_cols)
In [60]:
"""
---------------------------------------------------------------
Read in results from DEseq2 and format DFs nicer
---------------------------------------------------------------
"""
# Want to consider genes that are DE in FB & MB vs SC and HB and then look at those
# which are "consistently affected"
# i.e. for each gene make a "code of how it was affected"
df_all = pd.read_csv(f'{output_dir}df-all_epi-2500_{date}.csv')
df_dict = {}
deseq2_files = os.listdir(r_dir)
for f in deseq2_files:
if 'DEseq2' in f:
try:
de_df = pd.read_csv(os.path.join(r_dir, f))
de_df = de_df.rename(columns={de_df.columns[0]: 'u_id'})
gene_names = [s.split('-')[1] for s in list(de_df['u_id'].values)]
gene_ids = [s.split('-')[0] for s in list(de_df['u_id'].values)]
de_df['padj'] = de_df['padj'].fillna(1) # Replace Nan p values with 1.0s
# Replace Nans with 0's for other values
de_df = de_df.replace(np.nan, 0)
de_df[gene_id] = gene_ids
de_df[gene_name] = gene_names
df_dict[f] = de_df
except:
print(f)
Plot the venn diagrams of overlapping DE genes¶
In [61]:
df_dict.keys()
fb_unique_genes = []
mb_unique_genes = []
hb_unique_genes = []
sc_unique_genes = []
anterior_genes = []
posterior_genes = []
# DEseq2_CNS_wt_fb-hb_20210124.csv
fb_genes_all = []
mb_genes_all = []
hb_genes_all = []
sc_genes_all = []
cutoff = 2.0
neg_cutoff = -2.0
for k, d in df_dict.items():
if 'wt_fb-' in k:
gs = d[d['padj'] < 0.05]
gs = gs[gs['log2FoldChange'] > cutoff] # Only want those that are "selectively on " in FB
fb_genes_all += list(gs[gene_name].values)
if 'wt_' in k and 'mb' in k:
gs = d[d['padj'] < 0.05]
if k == 'DEseq2_CNS_wt_fb-mb_20210124.csv':
gs = gs[gs['log2FoldChange'] < neg_cutoff]
else:
gs = gs[gs['log2FoldChange'] > cutoff]
mb_genes_all += list(gs[gene_name].values)
if 'wt_' in k and 'hb' in k:
gs = d[d['padj'] < 0.05]
if k == 'DEseq2_CNS_wt_hb-sc_20210124.csv':
gs = gs[gs['log2FoldChange'] > cutoff]
else:
gs = gs[gs['log2FoldChange'] < neg_cutoff]
hb_genes_all += list(gs[gene_name].values)
if 'wt_' in k and 'sc' in k:
gs = d[d['padj'] < 0.05]
gs = gs[gs['log2FoldChange'] < neg_cutoff] # Only want those that are "selectively on " in FB
sc_genes_all += list(gs[gene_name].values)
In [62]:
# Now see which are also in the consistently affected genes
const_aff = pd.read_csv(f'{output_dir}df-consistent_epi-2500_{date}.csv')
const_aff_genes = const_aff[gene_name]
plot_venn([set(fb_genes_all),
set(mb_genes_all), set(hb_genes_all), set(sc_genes_all)],
['FB', 'MB', 'HB', 'SC'], 'SpecificGenes_loc', colours=[fb_colour, mb_colour, hb_colour, sc_colour])
plot_venn([set(fb_genes_all) & set(const_aff_genes),
set(mb_genes_all) & set(const_aff_genes),
set(hb_genes_all) & set(const_aff_genes),
set(sc_genes_all) & set(const_aff_genes)],
['FB', 'MB', 'HB', 'SC'], 'SpecificGenes_CAFF_loc', colours=[fb_colour, mb_colour, hb_colour, sc_colour])
Print out the genes for the table¶
This is just for supp. info
In [64]:
mb_caff = set(mb_genes_all) & set(const_aff_genes)
fb_caff = set(fb_genes_all) & set(const_aff_genes)
hb_caff = set(hb_genes_all) & set(const_aff_genes)
sc_caff = set(sc_genes_all) & set(const_aff_genes)
fb_spec = []
anterior_spec = []
for g in fb_caff:
if g not in mb_caff and g not in hb_caff and g not in sc_caff:
fb_spec.append(g)
if g in mb_caff and g not in hb_caff and g not in sc_caff:
anterior_spec.append(g)
mb_spec = []
for g in mb_caff:
if g not in fb_caff and g not in hb_caff and g not in sc_caff:
mb_spec.append(g)
posterior_spec = []
hb_spec = []
for g in hb_caff:
if g not in mb_caff and g not in fb_caff and g not in sc_caff:
hb_spec.append(g)
if g in sc_caff and g not in mb_caff and g not in fb_caff:
posterior_spec.append(g)
sc_spec = []
for g in sc_caff:
if g not in mb_caff and g not in hb_caff and g not in fb_caff:
sc_spec.append(g)
print(', '.join(fb_spec), '\n')
print(', '.join(mb_spec), '\n')
print(', '.join(hb_spec), '\n')
print(', '.join(sc_spec), '\n')
print(', '.join(anterior_spec), '\n')
print(', '.join(posterior_spec), '\n')
Run AP_axis functional¶
Run the functional enrichment and then re-read in the values from the GSEA and use this to plot the pathway terms.
In [19]:
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from adjustText import adjust_text
from sciviso import Vis
class Volcanoplot(Vis):
def __init__(self, df: pd.DataFrame, log_fc: str, p_val: str, label_column: str, title='',
xlabel='', ylabel='', invert=False, p_val_cutoff=0.05,
log_fc_cuttoff=2, label_big_sig=False, colours=None, offset=None,
text_colours=None, values_to_label=None, max_labels=10, values_colours=None,
figsize=(1.5, 1.5), title_font_size=8, label_font_size=6, title_font_weight=700):
super().__init__(df, figsize=figsize, title_font_size=title_font_size, label_font_size=label_font_size,
title_font_weight=title_font_weight)
super().__init__(df)
self.log_fc = log_fc
self.p_val = p_val
self.p_val_cutoff = p_val_cutoff
self.log_fc_cuttoff = log_fc_cuttoff
self.values_to_label = values_to_label
self.label_big_sig = label_big_sig
self.invert = invert
self.label_column = label_column
self.offset = offset
self.label = 'volcanoplot'
self.colours = {'ns_small-neg-logFC': 'lightgrey',
'ns_small-pos-logFC': 'lightgrey',
'ns_big-neg-logFC': 'grey',
'ns_big-pos-logFC': 'grey',
'sig_small-neg-logFC': 'lightgrey',
'sig_small-pos-logFC': 'lightgrey',
'sig_big-neg-logFC': '#df80ff',
'sig_big-pos-logFC': '#ffa366'} if colours is None else colours
self.xlabel = xlabel
self.ylabel = ylabel
self.title = title
self.figsize = figsize
self.max_labels = max_labels
self.values_colours = values_colours or {}
self.text_colours = text_colours or {}
def add_scatter_and_annotate(self, fig: plt, x_all: np.array, y_all: np.array,
colour: str, idxs: np.array, annotate=False, s=20):
x = x_all[idxs]
y = y_all[idxs]
ax = fig.scatter(x, y, c=colour, alpha=self.opacity, s=s, vmin=-10.0, vmax=10.0)
# Check if we want to annotate any of these with their gene IDs
if self.values_to_label is not None:
texts = []
labels = self.df[self.label_column].values[idxs]
for i, name in enumerate(labels):
if name in self.values_to_label:
lbl_bg = self.values_colours.get(name)
color = self.text_colours.get(name)
texts.append(fig.text(x[i], y[i], name, color=color, fontsize=6,
bbox=dict(fc="white", alpha=0.5, boxstyle='round,pad=0.1', lw=0)))
adjust_text(texts, force_text=2.0)
# Check if the user wants these labeled
if self.label_big_sig and annotate:
# If they do have a limit on the number of ones we show (i.e. we don't want 10000 gene names...)
max_values = -1 * self.max_labels
if len(y) < self.max_labels:
max_values = -1 * (len(y) - 1)
most_sig_idxs = np.argpartition(y, max_values)[max_values:]
labels = self.df[self.label_column].values[idxs][most_sig_idxs]
x = x[most_sig_idxs]
y = y[most_sig_idxs]
# We only label the ones with the max log fc
for i, name in enumerate(labels):
name = (' ').join(name.split('_')[1:]) #Format nicer
fig.annotate(name, (x[i], y[i]),
xytext=(0, 10),
textcoords='offset points', ha='center', va='bottom',
bbox=dict(boxstyle='round,pad=0.25',
fc=colour, alpha=0.2)
)
return ax
def plot(self):
"""
For annotation styling see: https://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.annotate
Returns
-------
"""
# if offset is not given, make the offset the smallest value in the dataset
if not self.offset:
vals = self.df[self.p_val].values
self.offset = np.min(vals[np.nonzero(vals)])
self.u.warn_p(['No offset was provided, setting offset to be smallest value recorded in dataset: ',
self.offset])
# x axis has log_fc, first only plot the values < cutoff
x = self.df[self.log_fc].values
y = -1 * np.log10(self.df[self.p_val].values + self.offset)
log_fc_np = self.df[self.log_fc].values
p_val_np = self.df[self.p_val].values
if self.invert:
x = -1 * np.log10(self.df[self.p_val].values + self.offset)
y = self.df[self.log_fc].values
sig_small_pos_logfc = np.where((p_val_np <= self.p_val_cutoff) & (np.abs(log_fc_np) < self.log_fc_cuttoff)
& (log_fc_np > 0))
sig_big_pos_logfc = np.where((p_val_np <= self.p_val_cutoff) & (np.abs(log_fc_np) >= self.log_fc_cuttoff)
& (log_fc_np > 0))
sig_small_neg_logfc = np.where((p_val_np <= self.p_val_cutoff) & (np.abs(log_fc_np) < self.log_fc_cuttoff)
& (log_fc_np <= 0))
sig_big_neg_logfc = np.where((p_val_np <= self.p_val_cutoff) & (np.abs(log_fc_np) >= self.log_fc_cuttoff)
& (log_fc_np <= 0))
# Plot the points
fig, ax = plt.subplots(figsize=self.figsize)
self.add_scatter_and_annotate(ax, x, y, self.colours['sig_small-pos-logFC'], sig_small_pos_logfc)
self.add_scatter_and_annotate(ax, x, y, self.colours['sig_big-pos-logFC'],
sig_big_pos_logfc, annotate=True, s=40)
# Negative
self.add_scatter_and_annotate(ax, x, y, self.colours['sig_small-neg-logFC'], sig_small_neg_logfc)
self.add_scatter_and_annotate(ax, x, y, self.colours['sig_big-neg-logFC'], sig_big_neg_logfc, annotate=True, s=40)
self.add_labels()
ax.tick_params(labelsize=6)
ax.tick_params(direction='out', length=2, width=0.5)
ax.spines['bottom'].set_linewidth(0.5)
ax.spines['top'].set_linewidth(0)
ax.spines['left'].set_linewidth(0.5)
ax.spines['right'].set_linewidth(0)
ax.tick_params(labelsize=6)
ax.tick_params(axis='x', which='major', pad=0)
ax.tick_params(axis='y', which='major', pad=0)
return ax
In [34]:
for f in os.listdir(fig_dir):
if 'GSEA' in f and 'KEGG' in f and 'padj' in f:
gsea_out = pd.read_csv(fig_dir + f)
gsea_out = gsea_out.fillna(1.0)
if len(gsea_out) > 1:
volcanoplot = Volcanoplot(gsea_out, 'NES', 'padj', 'pathway',
f.split("_")[2], 'NES', '-log10(p adj)',
p_val_cutoff=1.0, max_labels=5,
label_big_sig=True, log_fc_cuttoff=1.5, figsize=(1.5,1.5))
sns.set_style("ticks")
volcanoplot.plot()
save_fig(f'Volcano_gsea_{f.split("_")[2]}')
plt.show()
else:
print("Nothing significant for: ", f)
In [24]:
for f in os.listdir(fig_dir):
if 'GSEA' in f and 'Go' in f and 'svg' not in f:
gsea_out = pd.read_csv(fig_dir + f)
gsea_out = gsea_out.fillna(1.0)
volcanoplot = Volcanoplot(gsea_out, 'NES', 'padj', 'pathway',
f.split("_")[2], 'NES', '-log10(p adj)',
p_val_cutoff=1.0, max_labels=5,
label_big_sig=True, log_fc_cuttoff=1.5, figsize=(2,2))
sns.set_style("ticks")
volcanoplot.plot()
save_fig(f'Volcano_gsea_GO_{f.split("_")[1]}_{f.split("_")[2]}')
plt.show()
